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ABSTRACT 



A fully factorized two-dimensional Navier-Stokes flow solver has been 
developed and applied to the problem of predicting subsonic airfoil flutter in the 
light stall regime. The in viscid fluxes are evaluated with a central difference 
ADI scheme and fourth and second order numerical dissipation is used to obtain 
oscillation-free solutions. The performance of algebraic and one-equation 
turbulence models in predicting separated flow is explored for computing high 
Reynolds number steady flow and unsteady flows over an oscillating NACA 0012 
airfoil. Comparisons of the computed results with available experimental data 
indicate that even though the lift response is fairly well predicted, the 
computation of the pitching moment hysteresis loops is very sensitive to 
tuibulence modeling. Results computed with several current models are in good 
agreement whenever the steady stall angle is exceeded only slightly. However, 
they fail to capture the vortex shedding process leading to the onset of stall 
flutter. 
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I. INTRODUCTION 



A. BACKGROUND 

The development of numerical solution methods for two-dimensional 
Navier-Stokes equations during the past few years provides new tools for the 
investigation and prediction of airfoil flows. Of great interest is the study of flow 
separation on airfoils in unsteady motion which is usually referred to as the 
dynamic stall phenomenon. McCroskey et al [Ref. 1] performed a series of 
careful experiments which serve as valuable benchmark data. More recently, 
Lorber and Carta [Ref. 2] contributed important additional experimental 
dynamic stall information for a Sikorsky airfoil. They also investigated incipient 
torsional stall flutter [Ref. 3] and found that small-amplitude airfoil oscillations 
near static stall may be unstable. 

For a pitching airfoil the instantaneous work done on the fluid by the airfoil 
due to its motion is the product of the pitching moment about the axis of rotation 
and the differential change in angle of attack. This product usually is a positive 
quantity. However, if the net work per cycle of oscillation were to become 
negative then the fluid would be doing work on the airfoil. Once the airfoil 
begins to extract energy from the freestream, the amplitude of the oscillation 
will grow and finally diverge. This condition is known as stall flutter. Normally, 
flutter of an airfoil is due to a combination of torsion and bending. However, in 
this case flutter is caused by a single degree of freedom oscillatory motion. 
During a cycle of oscillation the lift coefficient and the pitching moment 
coefficient plotted versus angle of attack produce a hysteresis loop. It is the 
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pitching moment hysteresis loop that provides an indication of incipient stall 
flutter. As shown by Carta and Niebanck [Ref. 4], clockwise pitching moment 
loops represent negative aerodynamic damping and therefore cause oscillations 
of a free airfoil to grow' in amplitude, while counterclockwise loops cause such 
oscillations to decay. Hence torsional flutter will occur as soon as the area of 
the clockw'ise loop exceeds that of the counterclockwise loop. The 
aerodynamics of small-amplitude airfoil torsional oscillations near stall need to 
be investigated experimentally and computationally in order to examine the stall 
flutter mechanism. 

B. PURPOSE 

The first objective of this investigation is to test an unsteady, compressible 
Navier-Stokes code (N’sl.f) using an Alternating-Direction-Implicit (ADI) 
scheme based on the the Beam-Warming [Ref. 5] approximate factorization 
method to determine its ability to obtain realistic airfoil flow solutions for a 
variety of flow' regimes is examined. The accuracy of the numerical solution is 
investigated by comparing the computed solutions with available experimental 
data. Test cases include steady-state flow solutions at various flow speeds and 
angles of attack, as well as, unsteady flow' solutions over rapidly pitching and 
harmonically oscillating airfoils. In addition, the accuracy of the Xavier-Stokes 
solutions are further explored by comparing each case w'ith an unsteady, 
inviscid, incompressible, panel code (U2diif.f), and with a steady, 
incompressible, viscous/inviscid interaction code (Incompbl.O- The final and 
principal objective of this study is to determine the influence of mildly separated 
flow' during part of a harmonic oscillation cycle on blade stability. Results are 
presented for XACA 0012 airfoils showing the influence of various parameters. 
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Flow field details are also included in order to provide a better understanding of 
certain flow features that may lead to stall flutter. 
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II. GOVERNING EQUATIONS 



The continuity equation, the momentum equation, and the energy equation 
must be solved simultaneously in order to obtain a flow solution of a 
compressible viscous fluid about a body. A complete derivation of these 
equation can be found in various texts, eg., Anderson [Ref. 6]. However, main 
steps of the derivations are presented in the following sections. 

A. CONTINUITY EQUATION 

The continuity equation is the result of applying the physical principle of the 
conservation of mass to a finite control volume fixed in space. Simply stated 
the net flow out of a control volume through its bounding surface must be equal 
to the time rate of change of the mass inside the control volume. For any 
arbitrary control volume, the continuity equation can be expressed as 



B. MOMENTUM EQUATIONS 

The equation for the conservation of momentum is obtained by applying 
Newton’s second law, which state that the net force acting on a fluid particle is 
equal to the time rate of change of linear momentum of the fluid particle. In 




( 1 ) 



which for a two-dimension Cartesian Coordinate system becomes 



dp (9(pu) ^(pw) 

— — ' ’h _ 



d t dK dz. 



( 2 ) 
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Cartesian coordinates the momentum equations can be expressed as follows 



a 

and in the y-direction. 



d{pu) ^ (9(pu^ + p) ^ <9(puw) _ 

^ ^ ()x ^ 



(3) 



<9(pw) ^(puw 



Stress terms are as follows 



(4) 









2 (9w^ 

^ (9z y 



4xz=iU 






(5) 



C. ENERGY EQUATION 

The conservation law form of the energy equation is derived by applying 
the first law of thermodynamics (dE=dQ+dW) to a fluid particle. This leads to 



§ + ^|(E+pM + |[(E+p)w| = 

^ (pr„ + WT,, - Qx )+ ■£• (/^^XZ + WTxz - <ix ) 

where the total energy per unit volume and the heat flux are given by 



( 6 ) 



E = (e + l/2V ^)p, qi = -k dTjd^, 



(7) 



Here k is the thermal conductivity. 
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D. CONSERVATION FORM OF GOVERNING EQUATIONS 

The starting point for the numerical algorithm which is presented in the next 
section is the strong conservation law form of the two-dimensional Navier- 
Stokes equations. The non-dimensionalized vector form of the governing 
equations in conservation law form for a Cartesian coordinate system is; 



with, 



d t 9 k Re 



d¥ dG \ 
_v+ — y 

(5x ^ 



/ 



where. 



Tp 1 

'pu ' 



[pu 1 



fp 



w 



I pu p I I pwu 

<'=lpw^''=l. = 

Le J 



puw 
(e+p)u 



pw -l-p 



ro 1 



ro 1 



^ Fv=i;^ G.=| 



(e.p)w Ld LgJ 



( 8 ) 



(9) 



1 ^ 

■Wz 

y 



4 ( - 

= 3 - -w, 

'Tx^=P(Uz-wJ 

A ( 1 ^ 



fi = u T,, + WT„ -I-— ; -Ti 

4 " Pr(r-l) 



^ -a! 



Here the pressure is related to the conservation variables of q by the equation 



p = (7 ~ 0[e - 0.5p (u^ -I- v^)] 



( 11 ) 
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where y= 1.4 is the ratio of specific heats, and a is the local speed of sound 
given by, a^ = YP / p. 

The density and the velocities are non-dimensionalized with the freestream 
density poo and the free stream speed of sound aoo, respectively. The total 
energy is normalized by aoo^Poo- The time (t) scales as t* = t aoo / c, where (c) 

is a characteristic length, such as the chord length. The Euler equations are 
obtained from Equation (8) by dropping the viscous terms. The strong form is 
chosen because it enables shock capturing. 
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III. SOLUTION METHODS 



A. POTENTIAL FLOW METHOD (U2DIIF.F) 

For inviscid incompressible flow the conservation equations reduce to the 
Laplace equation which can be solved in a number of ways. We present here a 
brief outline of the widely used panel method for both steady and unsteady 
airfoil flow. 

1. Panel Method for Steady Airfoil Flows 

The frame of reference for the formulation of the steady flow problem 
is a fixed (x,y) coordinate system located at the leading edge of the airfoil. The 
freestream velocity is represented by +Voo and the x-axis passes from the 
leading edge through the trailing edge of the airfoil. 

In order to formulate a method for computing the flow around an 
airfoil, the airfoil surface is divided into (n) straight line segments or panels. The 
end points of the panels are called nodes, and since there are (n) panels there 
are then (n+1) nodes. Complex airfoil geometries can therefore be modeled 
with a greater number of nodes and panels. The trailing edge panel is the first 
panel, and the next panel continues on the lower surface in a clockwise manner 
until the n^h panel is reached at the trailing edge of the upper surface. 

Now that the frame of reference has been specified, two additional 
vectors must be defined, the unit normal vector ni, which is always 
perpendicular to the (i^^ ) panel and directed outward from the airfoil surface 
and the unit tangential vector ti, which is parallel to the (i^^) panel and is 
directed from the (n) node to the (n+1) node. 
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Two types of singularities are sufficient to model lifting airfoil flows. 
U2diif.f places uniform source distributions qj and a uniform vorticity distribution 
7 on each of the j-panels. Both singularities satisfy Laplace’s equation, a linear 
homogeneous second order partial differential equation. Since the solutions to 
linear PDE’s can be superimposed on each other, a simple flow can be added to 
another simple flow and so on, until a very complicated flow field is created. 

The flow field around an airfoil, represented by the velocity potential, 
can be constructed from the potential of the freestream flow added to the 
velocity potential of source and vorticity distributions, hence 

<^ = (p^+(p,+cp, 

where. 



= V„ +(x cosa + y sin a) 




summation of the three potentials yields 



(13) 

(14) 

(15) 



y»(x cosa+y sin a)+ [ rds- [ 0ds 

' ^ 2k J 2k ( 16 ) 

Equation (16) can now be evaluated at any point in the flow field by evaluating 
the integrals along the airfoil contour (s), where the flow field point is located at 
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a distance (r) at an angle (0) measured from a point on the airfoil. Introducing 
n panels and then summing the contributions of each panel yields 

d>=V„ + (x cosa + y sina) + 

Once the boundary conditions are defined the system of (n) equations can be 
solved for (n+1) unknowns. 

Next, it is useful to introduce the concept of influence coefficients. An 
influence coefficient is the velocity induced at a point in the flow field (field 
point) by a unit strength singularity, source or vorticy, placed anywhere within 
this field. U2diif.f places these points on each panel. A detailed description of 
the use of influence coefficients is found in Teng [Ref. 7]. 

For the steady flow problem the influence coefficients are given in 

Table 1 . 

TABLE 1. INFLUENCE COEFFICIENTS 



AHij 


Normal component induced at i^*^ control point by 
unit source distribution on panel. 


Atij 


Tangential component induced at i^h control point 
by unit source distribution on panel. 


BHij 


Normal component induced at i^^ control point by 
unit vorticity distribution on panel. 


B>ij 


Tangential component induced at i^h control point 
by unit vorticity distribution on jth panel. 


6ij 


Angle of control point (i) and panel (j) measured 
from x-axis. 
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There are only two boundary conditions that must be satisfied in order 
to solve the lifting airfoil flow problem. The first boundary condition is the 
condition that the flow remains tangent to the airfoil surface and the second 
condition is that the pressures on the upper and lower trailing edge panels must 
be equal. This condition is known as the Kutta condition. Pressure is related to 
velocity through Bernoulli’s equation for steady potential flow, which allows the 
Kutta condition to be expressed as 



' 'lower ' 'upper 



The flow tangency condition is expressed as 



(18) 



(v"°™*'). =Q i = l,2,.., n 



(19) 



Equation (18), the Kutta condition, using the influence coefficient concept is 
expressed as follows 



-X[A|jqJ-yZ[B;j|-V„cos(«-e,) = 

j=l j=l 

i(A;,qJ+ ri|B'J+ v.cos(a - ej 

j-1 j-1 (20) 

Equation (19) the flow tangency condition, becomes, using the influence 
coefficient form 



ik"qJ+rS[B;J + V.sin(a-eJ=a i=l,2,...,n 

j=i j=i (21) 
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These equations can then be written in the following matrix form 




( 22 ) 



and solved with the method of Gauss Elimination with Partial Pivoting. 

2. Unsteady Numerical Formulation 

A brief and concise description of the unsteady numerical formulation 
is found in Krainer [Ref. 8]. To model unsteady flow around an airfoil using the 
panel method, N unknown source strengths and one unknown vorticity strength 
are located along the airfoil, a wake panel of unknown vorticity strength, length 
and orientation is attached to the trailing edge of the airfoil. This makes a total 
of N+4 unknowns. The flow tangency requirement at each of the N panels 
represents a system of N equations. Using influence coefficients and using the 
subscript k to count time, the tangential flow condition of the i^^ element is 
written as follows 



b;. +[(V. + (U(t)i + V(t)j)+ Q(i)(xi - >])) n.] 

j-1 j=l k 

+(rwUBr,.,l + I|(c:,)(r^, - r„) I =0 i = 1,2 n 

rn=l (23) 

where Voo is the mean velocity, (U(t)i + V(t)j) is a time dependent translational 
velocity and Q(t) is a rotational velocity. The vectors i and j are in the airfoil 
fixed coordinate system. 
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The Kutta condition is a single equation that determines the circulation 
around the airfoil. In this case the pressures at the upper and lower surface are 
equated at the midpoints of panels adjacent to the trailing edge (i.e., the first and 
last panel) 



[(vl)J +[(vn)J -2j 



2 JdT^ 






= 2 



7k -7k-i 

*^k ^k-i ^ 24 ) 



The Helmholtz theorem is another single equation that provides a 
relation for the strength of the vorticity distribution along the wake element. 
The Helmholtz theorem states that any change in circulation around an airfoil 
must be countered by a change in vorticity in the wake of equal magnitude but 
of opposite sign. This theorem can be written as follows 



Ak(7w)k = ^(7k-i-7k) 



(25) 



where / is the length of each individual wake element and y its vorticity 
strength. Therefore, the vorticity in the wake element is equal to the negative 
change in circulation around the airfoil with respect to the k-1 timestep. 

To acquire the two additional equations required to solve for the N+4 
unknowns, assumptions about the geometry of the wake panel are made. First, 
the wake panel is oriented in the direction of the local resultant velocity at its 
midpoint, as viewed in the (moving) airfoil frame of reference. 




( 26 ) 
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(vyw)k and (v’^w)k are the x- and y-velocity components at the midpoint of 
each wake panel. Second, the length of the wake panel is assumed 
proportional to the magnitude of the local resultant velocity at its midpoint and 
to the timestep. 

Ak = ^(v:)J +[(v:)J (tk -tk_i) 

(27) 

The nonlinearities in the Kutta condition and in the wake panel 
assumptions necessitate an iterative solution procedure. 

B. VISCOUS/INVISCID INTERACTION METHOD 

(INCOMPBL.F) 

The theory and numerical methods presented in this section are taken from 
the work of Cebeci and Bradshaw [Ref. 9] and the investigations performed by 
Krainer [Ref. 8] and Snir [Ref. 10] As stated in the introduction the results 
from the Navier-Stokes code (Ns2.f) were compared with a viscous/inviscid 
interaction method code (Incompbl.f) made available by Cebeci. An 
abbreviated description of this method is given in this section. For a complete 
description see the original publications, References 8 and 9. 

The governing principle behind the viscous/inviscid interaction method is an 
approximation to the Navier-Stokes equation that allows a flowfield to be 
divided into an inner viscous region and an outer inviscid region when the 
Reynolds number is sufficiently large. This concept of a thin viscous layer near 
the surface of a body and an outer region where these viscous effects are small 
compared to the inertial effects is known as the boundary layer theory. Unlike 
the Laplace equation, which is the governing equation for potential flow, the 
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thin shear layer equations are nonlinear 



d\x dw 
dx dy 



u 



dx 





with the boundary conditions 



( 28 ) 



y=0, u=v=0, 

y=oo, u=Ue(x) (29) 

The direct boundary layer method solution to the boundary layer equations 
by the Keller Box Method involves four steps. First the boundary layer 
equations are transformed into a system of first order differential equations. 
Next the Keller box method is used to approximate the first order differential 
equations by simple centered differences and two-point averages, using values 
at the corners of one difference molecule only. Newton’s method is used to 
linearize the resulting algebraic equation. The Keller block elimination method 
is then used to solve the resulting block tridiagonal system. The procedure 
described above does not provide solutions to flows that are separated or have 
regions of reverse flow. 

1. Interaction Method 

The interactive boundary layer method provides a special coupling 
between the inner viscous flow and the outer inviscid flow, which enables 
reverse and separated flows to be calculated. In such areas, the external 
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velocity is substantially changed by the viscous effects and can no longer be 
considered as a known boundary condition for the boundary layer flow. 

The general approach to the solution is the same as for the direct 
method with modifications. Since the outer flow is unknown, the velocity at the 
edge of the boundary layer is given by the interaction law and is written as 



where Ue(x,ye) is the total velocity at the edge of the boundary layer, Uei(x) is 
the velocity computed by the inviscid panel method, 5* is the displacement 
thickness, and the integral term is known as the Hilbert integral. 

The following transformations are used to transform the boundary 
layer equation into a system of first order differential equations 




(30) 




(31) 



The boundary layer equation takes the form 



f=U, U'=V, W’=0, 




(32) 



with boundary conditions: 



ri=0. 



U(x,0)=0, f(x,0)=0, 

U(x,Tle)=W(x,Tle) 






(33) 
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The velocity at the edge of the boundary layer now becomes 




The finite difference box method is used to solve the equations, in the 
same way it was used for the direct case, but with two additions. First in areas 
of flow reversal the u3u/3x is omitted to assure stable integration. And second, 
the edge velocity is approximated by the relation presented in the next section. 

By using central differencing to approximate the differential equations, 
a system of nonlinear algebraic equations is obtained for the unknown variables 
(which are f, U, V, and W). To solve the system of equations, the system is 
linearized by the Newton iterative procedure, and the resulting linear system is 
solved for the new unknown variables which are the increments 5f, 5U, 5V, and 
5W. 

The solution procedure is repeated until the change in the increments 
is negligible compared to the preceding iteration. The iterative process is 
performed again at the next downstream station. 

2. Interactive Model 

The interactive model is used to couple the boundary layer to the 
external flow. It is needed in areas where strong interaction occurs, and both 
the boundary layer and the outer flow must be solved simultaneously. The 
external velocity is assumed to consist of a potential flow term and a correction 
term due to viscous effects. [ue(x)=uel(x)+ue5(x)j 
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The normal velocities at the surface of the airfoil, induced by sources 
distributed on the surface, displace the streamlines from the surface in the same 
way that the actual boundary layer displaces them. 

Several assumptions are made in order to express the correction term 
in the form of the Hilbert integral. First, the surface is approximated to be a flat 
plate, and the normal velocity is then half the local source strength cr(x). 
Second, the inviscid velocity does not change across the boundary layer’. 
Therefore, the local horizontal velocity induced by the source distribution, is the 
correction term to the inviscid velocity, and can be represented by the Hilbert 
integral 



The integration is carried out on all the sources on the surface. The Hilbert 
integral is then approximated by a finite series 



where cik is a matrix of interaction coefficients which are functions of the 
geometry only. 

Since the computation of Ue6 involves values of 5* downstream of the 
current x location, which are not known yet, these terms are taken from the 
previous iteration using a relaxation formula. 

3. Turbulence Model 

The turbulence model used in the code is the Cebeci-Smith model 
[Ref. 11] and Michel’s method is used to predict the transition from laminar to 
turbulent flow. 




(35) 




(36) 
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C. EULER AND NAVIER-STOKES CODE (NS2.F) 



In order to facilitate the numerical implementation for arbitrary complex 
flow domains, the Navier-Stokes equations are transformed to a generalized 
coordinate system, (^,0, using the transformations 

=(^(x,z) 

C = C(x,z) (37) 

The grid spacing in the curvilinear space is uniform and of unit length so that 
unweighted differencing schemes can be employed for the numerical 
implementation. The grid points in the Cartesian system, referred to as the 
physical domain, correspond via a one to one relationship to the points of the 
curvilinear transformed system, referred to as the computational domain. 
Singularities of the transformation may occur on the computational boundaries. 
Such singularities occur for grids with multiple connected regions. The 
transformation of the governing equations from the physical domain to the 
computational domain is obtained in most cases by numerically evaluating the 
metrics and the Jacobian of the transformation. The derivatives with respect to 
the X, y variables can be expressed in terms of the new variables by the chain 
rule. 



d ^ d y d 
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The metrics and the Jacobian of the transformation are given by 



= J z;, Cx = -J z^, 

~ "J Cx = J X^, 

(x|Z{-XjzJ 

The Navier-Stokes equations written in generalized curvilinear coordinates 
retain the strong conservation law form expressed by Equation (8). The 
governing equation for generalized curvilinear coordinates is 

^ A 

The conservation variable vector <1 and the inviscid fluxes F and G in the 
transformed space are given by 

\p 1 Tpu 1 Tpw 1 

._J_lpul + I ^_l_|pWu+C,p I 

^ jipwl’ jjpwU+4p j jjpwW + ^^p j 

[e J |_(e+p)U-4pj + (43) 

where U and W are the contravariant velocity components along the ^ and 
directions respectively, given by: 



(39) 

(40) 

(41) 



U = + W = Ci + C,u+CzW 

The viscous flux terms are transformed as 

Fv = J (4Fv + 4G v), G , = y (C.F, + CG, ) 



(44) 



(45) 
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with the stress terms from Equation (10) expressed in terms of the transformed 
variables ^ as 



'Z’xz = + CzU J+ + C,w J] 

^zz =3A^[(4w4+CzwJ-^fe,u^+Cxujj 



Ai 






^ dsL^ ^ d?i 



fl r (5a^ 



(46) 



When the thin layer approximation [Ref. 12] to the two-dimensional 
conservation law form of the governing equations is applied, they take the 
following form 



^ ^ dG 1 1 

d ~Rt[dC \ 



where, 



(47) 



rp 1 


rpu 1 


llpu 1 


^_j_|puu + 4p 


J 1 pw 1’ 


J j pwu 4 p I 


[e J 


L(e4-p)U-^pJ 



fpw 1 

. ] I pWu + I 

^-jjpwW-h^.p j 

L(e + p)W-CpJ 



(48) 
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The viscous flux term is transformed as 



^ /X m,u^ + (^/3)m2 C, I 

''~J ;x mjW^ + (/x/3)m2 

m,m3 + (;x/3)m2m4^ 



( 49 ) 



Here, 



m2=CxU^+CzW^ 

m4 = CxU + CzW 



1 I" (9a^ 1 



(50) 



1. Numerical Grids 

In order to compute flow solutions of partial differential equations 
(such as the Navier-Stokes Equations) with finite differences, a discretized 
version of the physical domain must be generated. During the course of this 
investigation, several different NACA 0012 airfoil C-type grids were generated 
using two different grid generation methods. A grid generation code call 
GRAPE that uses the Poisson differential equation was used in the early stages 
of investigation. The inviscid and viscous grids shown in Figures 1 and 2 were 
generated from an algebraic grid generation code. Note the finer resolution of 
the viscous grid near the airfoil surface. This resolution is required to 
accurately represent the boundary layer. The flow solutions presented in 
Chapter IV were computed by using the grids shown in Figures 1 and 2. The 
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advantage to using an algebraic grid , vice a PDE method or a conformal 
mapping method, is that it is computationally efficient. 

A simple description of the algebraic grid generation method is offered. 
First the airfoil surface coordinates are unwrapped to form a simple curve, by 
separating the lower trailing edge from the upper trailing edge. This curve is 
the starting point for the generation of the computational domain. Lines are 
drawn normal to the curve and grid points are then located at intervals on these 
lines. Lines may be clustered near the leading or trailing edge of the airfoil; or 
at some other area of interest, eg., a shock location. The resolution of the 
points normal to the surface is achieved by the use of an algebraic function that 
provides uniform stretching normal to the body surface. Once the grid is 
created in the computation domain it is then transformed back into the physical 
domain using inverse transformations. 
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Figure 1. NACA 0012 Inviscid 151x64 C-Grid 
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Figure 2. NACA 0012 Viscous 161x64 C-Grid 
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2. The Numerical Scheme 

The integration in time is performed implicitly with a second order 
accurate scheme of the form 



^ ■‘i )+T— +0 



2 di. 



j 



(51) 



Using this trapezoidal rule differencing scheme to approximate the time 
derivative of the conservative dependent variable q vector, the unknown 
conservative variables at the n+1 timestep are given by 



. n n + 1 , 

Aq =q +y5 



' dt 









(52) 



The spatial derivatives of the governing equations are discretized using central 
differences, and the right hand side of Equation (52) becomes 



Aq" = -y{(5jF + 5jd )" + (SfF + Sfi )"' } 



2 1" * - " > ' j (53) 

The nonlinear terms F and G at the n+1 timestep are linearized as follows 



F"^' =F"+_Aq"^‘ +0(At^) = F" + A"Aq""‘ +0(At") 

^ (54) 

where A" is the flux Jacobian matrix given in Appendix A. Substitution of the 

linearized form of the flux vectors in Equation (54) yields 
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The two-dimensional operator on the left hand side of Equation (55) is 
approximately factorized as follows: 



Lll + f + [Aq-; = At{5,F",-5,G:,} 



2 ^ 



(56) 



On the right hand side of Equation (56) an explicit dissipation term Dexpl, is 
added for numerical stability. In order to obtain an oscillation free solution a 
second order, dissipation term Dimpl is added to the left hand side spatial 
operators. Thus, the complete discretized form of the governing equations 
becomes 



{ni+f(^iA:,.+(D,jJ} 




(57) 




>:,J = (RHS)" 


(58) 



This equation is solved by performing two sweeps as follows: 



III + (d„ J lUq'u' = (RHS)" 



(59) 






(60) 



During each sweep the following linear block tridiagonal system of equations is 
solved 



+b..,Aq,,, 



(61) 
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where Aq = Aq*^^+1 and f = (RHS)*^ when the ^ direction sweep is performed, 
and Aq = Aq^^+l and f = (RHS)*fi when the ^ direction sweep is performed. 
The block matrices have the form 



At ^ . ^i.k 

^i.k “ 9 ^i.k t 

^ -^i+Uk 


(62) 




(63) 


At ^ J i it 




Ci.k = ^ ■ 

^ i-1 ,k 


(64) 


Boundary Conditions 





All flows were computed at subsonic free stream speeds. For 
subsonic inflow and outflow boundaries the flow variables are evaluated using 
zero order Riemann invariant extrapolation. At the inflow boundary there is one 
incoming and three outgoing characteristics, therefore, three variables, density 
(p), normal velocity (wi), and pressure (p) are specified and the fourth variable 
axial velocity (ui), is extrapolated from the interior. The inflow boundary 
conditions are given by 



(65) 



|r-.| 
Pi = 

bij 



w, =w^ 



S, 



Pi = 



r y 
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where R"^!, R'2 are the incoming and outgoing Riemann invariants given by 
R;=u. + 2a./(r-l). R- = u,-2aj/()--l) 

( 66 ) 

At the outflow boundary there is one incoming and three outgoing 
characteristics and only one quantity, pressure, is specified while the others are 
extrapolated from the interior. For the density and the normal velocity, simple 
first-order extrapolation is used, and the axial outflow velocity is obtained from 
the zero order outgoing Riemann invariant. The outflow boundary conditions 
are given by 



A =P2 

u, = R; -2a,/(r- 1), a, = VWA 

Wj = W2 

Pi ~ P 2 

(67) 

On the body surface the nonslip condition is applied for the velocities. The 
density and pressure are obtained from the interior by extrapolation. For the C- 
type grids used in this study averaging of the flow variable at the wake cut is 
used. 

4. Turbulence Models 

a. Baldwin-Lomax Turbulence Model 

The Baldwin-Lomax model [Ref. 12] is a two -layer, inner and 
outer eddy viscosity model for the computation of two- and three-dimensional 
flows. Patterned after the Cebeci-Smith model [Ref. 11], it incorporates a 
modification that bypasses the need for finding the edge of the boundary layer. 
This is achieved by introducing the vorticity in place of the boundary layer 
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thickness. The turbulent effect is simulated by an inner eddy viscosity, given 
by: 






y ~y crossover 



(68) 



where ycrossover is the smallest value of y at which the values of the inner and 
outer eddy viscosities are the same and the Prandtl mixing length 1 is 



-y^l 

jj 



r 

l = kyi 1 —exp 



(69) 



The magnitude of the vorticity in two-dimensions lo)l and y+ are defined as 
follows 



N= 



^ <9v'] 



_ PwUry 



Jp^ 



and ' Pw Pw ( 70 ) 

where A+ is a constant. The subscript w denotes values at the wall or airfoil 
surface in this case. 

The outer eddy viscosity is given by the following expression 



(Pl)outer= '9:ciPF(y)wAK^(y) 



KLEB 



' crossover 



^y 



(71) 



where k and Ccp are constants and F(y)wake=ymax Fmax for boundary layers and 
F(y)wake=Cwk ymax (UDIF^/Fmax ) for wakes and separated boundary layers. 
The quantities ymax and Fmax are determined from the expression 



r ( 

F(y) = y \<4 1^1 -exp 



A* JJ 



(72) 



30 



The exponential term of Equation (72) is set to zero for wakes. Fmax the 
maximum value of F(y) that occurs in a profile and y^ax is the value of y at 
which it occurs. F(y)KLEB is the Klebanoff intermittency factor given by 



The quantity Udif is the difference between Umax=Uy=ymax and the minimum 
total velocity in the profile 



The second term in Udif is taken to by zero except in the case of wakes. 

The constants were determined to achieve agreement with 
Reference 11 and can be found in the original document. This model 
completely models turbulent flow over an airfoil, it is relatively easy to code, it 
does not degrade the solution convergence, and it is computationally efficient. 
b. Johnson-King Turbulence Model 

The Johnson-King model [Ref. 13] takes into account the 
convective and diffusive effects on the Reynolds shear stress -u’w’ in the 
streamwise direction. The eddy viscosity is given by 



where Vt , Vto describe the eddy viscosity variation in the inner and outer part of 
the boundary layer. The inner eddy viscosity is computed as 




(73) 




(74) 




(75) 
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2 _fy/y^+ j 

max , where D =1 -e , where the constant A+=15. The 



outer viscosity is given by 



Vto=o(x)[0.168Ut5xy] 



(76) 



where yis the Klebanoff intermittency function y=[l+5.5(y/5)6]‘l and a(x) is the 



-u’w’lmax along the path of maximum shear stress. This model accounts for the 
effects of convection and diffusion on the Reynolds stress development through 
the solution of the following ODE 




here Cdif and ai are modeling constants, Um is the maximum average mean 
velocity and g=[-u’w’lmax]'^'^^, geq=[-u’w’lmax,eg]'^^^ where Lm is the 
dissipation length evaluated as 



solution of the ordinary differential equation which describes the development of 



L„ = Q40y >v/5<Q225 

L„ = Q09y y„/5>Q225 



(78) 
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The equilibrium shear stress geq in Equation (77) is determined from the 
following equilibrium eddy viscosity distribution 




=DV/-u'w')„^ 
''.0.1 =<J(x)(Q168U,5>| 



(79) 



An implicit Euler method is used for the numerical solution of Equation (77), 
and the maximum shear stress at each iteration level is updated as follows 



Solutions with the Johnson-King turbulence model were obtained as follows. 
First a convergent solution using the Baldwin-Lomax turbulence model for the 
entire flow field was computed. Then the Johnson-King model was applied only 
to the upper part of the airfoil. To initiate the solution a(x) in Equation (76) is 
set equal to unity and it is allowed to change according to Equation (80) until 
the final solution is obtained. 

c. Algebraic RNG-based Turbulence Model 



e model based on Renormalization Group (RNG) Theory of turbulence [Ref. 14] 
were proposed for the closure of the Reynolds-averaged Navier-Stokes 
equations. The algebraic model, although free from uncertainties related to the 
experimental determination of empirical modeling constants, still requires 
specifications at an integral length-scale of turbulence which reduces the 



Cj(x)""‘ = ct(x) 




(80) 



Recently an algebraic eddy viscosity, as well as a two-equation k- 
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generalities of the model. Here the integral scale is assumed proportional to the 
distance from the wall and the eddy viscosity for the RNG-based algebraic 
turbulence model is obtained as in Martinelli [Ref. 15] from the following 
formula 



1 




where v=Vt+vi, H is the Heavyside step function and ([) is the dissipation 
function (})=Xij ( du[ / Oxj ). The RNG-based turbulence model is applied only for 
the suction surface separated flow region while the pressure side and the wake 
regions are computed with the Baldwin-Lomax model. 
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IV. RESULTS AND DISCUSSION 



A. METHOD OF INVESTIGATION 

As stated in the introduction this investigation is divided into three parts; 
validation of the Navier-Stokes code (Ns2.f) through comparison with 
experimental data and a well tested inviscid unsteady Panel method code as 
well as a steady viscous/inviscid code. Finally evaluation of the effects of 
reduced frequency, mean angle, amplitude, Mach number, and Reynolds 
number on unsteady flow behavior are examined. Also, during the course of 
this study it became necessary to investigate the effects of several turbulence 
models on the results computed for the harmonically oscillating cases. 
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The following table lists the cases studied during the course of this 
investigation; 



TABLE 2. TEST CASES 



Case 


Type of Motion 


Parameters 


1 


Steady State 
(Transonic) 


M=0.7, a=1.4°, Re=9xl06 


2 


Steady State 
(Transonic) 


M=0.799, a=2.26°, Re=9xl06 


3 


Steady State 


M=0.3, Re=3xl06 


4 


Steady State 


M=0.3, Re=4xl06 


5 


Rapidly Pitching Airfoil 
(Ramp) 


M=0.3, k=0.01272, Re=2.7xl06 


6 


Harmonically Oscillating 
(Multiple Parameters) 


a(t)= Ao°±A 1 °sin(o)t) 


7 


Harmonically Oscillating 
(Small Amplitude) 


a(t)=13°±2.5°sin(o)t), k=0.2, Re=4xl06 


8 


Harmonically Oscillating 
(Small Amplitude) 


a(t)=13“±2.5“sin(cot), k=0.1, Re=2xl06 


9 


Harmonically Oscillating 


a(t)=9°±5°sin((ot), k=0.2, Re=4xl06 


10 


Harmonically Oscillating 


a(t)=10±5.5“sin(o)t), k=0.1, Re=4xl06 


11 


Harmonically Oscillating 


a(t)=ir±5“sin((ot), k=0.1, Re=4xl06 
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1. Steady-State 

Flow solutions can be obtained either by rotating the grid or by rotating 
the flow. Cases 1 and 2 solutions were obtained by rotating the grid. 

First the procedure was used for the computation of steady flow over 
the NACA 0012 airfoil. Steady-state solutions were obtained for subsonic and 
transonic flow regime. Flow solutions for Case 1 and 2 (Table 2) were 
computed by rotating the direction of the flow relative to the stationary grid. 
However, in Cases 3 and 4 the grid was rotated to the desired angle using a 
simple Fortran code called rotgr.f and the free stream flow remained parallel to 
the x-axis. The program also translates the grid to a desired pivot point. In this 
thesis the pivot point was always chosen to be the quarter chord point. Either 
method produced the same solution. The steady-state solutions tested the 
accuracy of Ns2.f and provided initial solutions for the unsteady cases. Ns2.f 
was also compared to an existing steady, incompressible, viscous/inviscid 
method (Incompbl.f) obtained from Dr. T. Cebeci, California State University at 
Long Beach. 

a. Case 1. M=0.7, a=1.4'‘, Re-9xl0^ 

The flow conditions for the first test case are Moo=0.7, at an angle 
of attack of 1.4 degrees, and a Reynolds number of nine million. Explicit 
boundary conditions were used and the Baldwin-Lomax turbulence model was 
selected. The solution was obtained on a 161x64 point grid. This case was 
chosen because the solution converged rapidly. Figure 3 shows good 
agreement between the inviscid (Euler) pressure distribution and the measured 
values of Harris [Ref. 16]. Next, a Navier-Stokes solution was obtained whose 
convergence history is given in Figure 4. The solution was carried out to 2500 
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iteration at a timestep of approximately 0.08. The pressure distribution 
computed by Ns2.f (Figure 5) is also in good agreement with the measured 
values. However, the suction peak is more accurately predicted by the inviscid 
solution. For a free stream Mach number (Moo) of 0.7 a maximum Mach 
number of 1.1 (Figure 6) is obtained on the upper surface near the leading edge. 
The density contour plot, Figure 7, shows the smooth contours indicative of a 
fully converged solution. The velocity field. Figure 8, shows the boundary layer 
thickness, represented by the thick dark line. The sonic line is represented by 
the line extending from 7% chord to 20% chord. 
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Cp plot for 00 1 2 M=0.7 a= 1 .4 
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Figure 3. Pressure Distribution (Inviscid), M=0.7, a=1.4°, Re=9xl0^ 
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Chordwise Location, x/c 



Convergence History, M=0.7, a=1.4, Re=9x 10**6 
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Figure 4. Convergence History, M=0.7, a-1.4 , Re-9xl0^ 
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Iteration 



M=0.7, a= 1.4 deg, Re=9x 10**6 




Figure 5. Pressure Distribution , M=0.7, a=1.4°, Re=9xl0^ 
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Chordwise Location, x/c 
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Figure 6. Mach Contour, M=0.7, a=1.4°, Re=3xl06 



42 



DENSITY 



CONTOUR LEVELS 
0.62000 
0.64000 
0.66000 
0.68000 
l ).’/(/000 



0.700 
1.40 DEG 
9.00x10*^6 
16. 

16k64 






i,'.S8*KX') 




1 .* )8XH)n 
1.10000 
1.12000 
1.14000 
1.16000 
1.18000 
1.20000 
1 .22000 
1.24000 
1.260un 



MACH 

ALPHA 

Re 

TIME 

GRID 



Figure 7. Density Contour, M=0.7, a=l,4°, Re= 9 xl 06 



contour leuels 

H. RROfin 
U . bOOUO 
It . synni ) 

II. '.-II Mil 

i * f.i ;ni 



0. 700 


9006 


1 . ^«EG 


8LPHR 


9>OOjcJD 


• *6Re 


IG. 


nnc 


6£4 


GRID 



I . 'I -lii 

.) h-JIlUll 
11 . I.UTII) 
;j. biu iin 

il. /Ill Uil 
I, 

'LUl 



1 . .-; 4 rnn 
0,86000 
0.00000 
0. 90000 
0.92000 
u. y-iuuu 

0. 96000 

0. 9H0D0 

1 . uuooo 




Figure 8. Velocity Field, M=0.7, a=1.4°, Re= 9 xl 06 
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b. Case 2. M=0.7 99, a=2.26\ Re=9xl0^ 

The next steady-state case investigated was at a Mach number 
equal to 0.799, at an angle of attack of 2.26 degrees and a Reynolds number of 
nine million. As in Case 1, an Euler solution was first obtained and compared to 
measured data as well as to a viscous solution obtained with an upwind scheme 
using the Johnson-King turbulence model. As seen in Figure 9, the pressure 
distribution for the inviscid solution is in fair agreement with the measured data 
until the shock is reached at the 50% chord point. After this point the computed 
solution fails to predict the pressure jump due to the presence of a shock at 
mid-chord. Note that the upwind scheme in combination with the Johnson-King 
turbulence model accurately predicts the location of the shock, even though the 
magnitude of the pressure gradient is off slightly. Next, a viscous solution was 
computed by Ns2.f with the Baldwin-Lomax turbulence model. Fourth order 
dissipation was added to produce a smooth flow solution across the weak shock 
located at mid chord. The solution was continued for 6000 iterations. The 
spike in the convergence history (Figure 10) was due to the addition of more 
fourth order dissipation. It is seen that the shock location was not predicted. 
At this angle of attack and Mach number the flow starts to separate and the 
Baldwin-Lomax turbulence model fails to model this separated flow region. 
Figure 1 1 shows the viscous pressure distribution comparison. The sonic line is 
shown on the Mach contour plot of Figure 12. A maximum Mach Number of 
1.45 is reached in this supersonic flow region. Note the steep Mach gradient 
(Figure 12) at approximately 60% chord. Figures 13 and 14 show the pressure 
contour and density contour, respectively. Again the steep density gradient 
indicates the shock location at approximately 60% chord vice the measured 
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location at 50% chord. Finally, the velocity field and sonic line are shown in 
Figure 15. 



NACA0012 M=0.799 at a=2.3 de 
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Figure 9. Pressure Distribution (Inviscid) , M=0.799, a=2.26% 

Re=9xl06 
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Chordwise Location, x/c 



Convergence History, M=0.799. a=2.26, Re=9x 10**6 
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Figure 10. Convergence History (Ns2.f)? M=0,799, a=2.26°, 

Re=9xl0^ 
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Iteration 



M=7.99, a=2.26 deg., Re=9x 10**6 
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Figure 11. Pressure Distribution, M=0.799, a=2.26°, Re=9xlQ6 
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Chordwise Location, x/c 
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Figure 12. Mach Contour, M=0.799, a=2.26°, Re= 9 xl 06 
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Figure 13. 



Pressure Contour, M=0.799, a=2.26% 



Re=9xl06 
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Figure 14. Density Contour, M=0.799, a=2.26°, Re= 9 xl 06 
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Figure 15. Velocity Field, M=0.799, a=2.26% Re=9xl06 
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c. Case 3. M=0.3, Re=3xl0^ 

In order to produce C] vs a and Cm vs a curves, solutions were 
computed for a=2°, a=4°, a=6°, a=8°, a=9°, a=10°, a=12°, a=13°, a=13.5°, 
a=14°, and a=15°. However, only results for a=4°, a=ll°, and a=14°, are 
presented. 

Figure 16 shows the convergence history of Ns2.f for a 161x64 
algebraic C-type grid. The Courant-Friedrichs-Levy (CFL) condition was 
approximately 2000-2500 which corresponded to a timestep between 0.005 and 
0.01. A comparison of the computed pressure distributions at a=4° (Figure 17) 
shows good agreement between the Navier-Stokes calculations using the 
Baldwin-Lomax turbulence model and the viscous/inviscid code predictions. 
The measured suction peak is somewhat higher than predicted by the codes. 
Figures 18 and 19 show fair skin friction agreement. Note that the Navier- 
Stokes code assumes turbulent flow over the entire airfoil, while the 
viscous/inviscid interaction code uses Michel’s method to predict the 
laminar/turbulent transition. As seen in Figure 18, the transition point for the 
upper surface is at 10% chord. A Mach contour plot is shown in Figure 20 and 
a density contour plot with the normalized stagnation pressure (0.98) contour 
overlayed in Figure 21. A maximum Mach number of 0.48 was attained. The 
density contours are smooth with steep gradients at the leading edge. 

In Figures 22 through 29 similar results are presented for angles of 
attack of 11 and 14 degrees. It is seen that the agreement between the two 
codes and the measured pressure distributions is quite satisfactory, although the 
suction peak tends to be underpredicted by both codes. 



54 



Figures 30 and 31 show the Ci vs a and Cm vs a curves, 
respectively. Stall is seen to occur at approximately a=13.5°. The Navier- 
Stokes code reproduces the experimental lift and moment values up to the 
maximum lift value quite well, whereas the viscous/inviscid code shows greater 
deviation for angles of attack exceeding about eight degrees. 
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Convergence History, M=0.3, a=4.0, Re=3x 10**6 







Figure 16. Convergence History, M=0.3, a=4°, Re= 3 xl 06 
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Iteration 



Steady State: M=0.3, Re=3x 10**6 
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Figure 17. Pressure Distribution, M=0.3, a=4°, Re=3xl0^ 
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Chordwise Location, x/c 



M=0.3, a=4.0 deg, Re=3x 10**6 
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Figure 18. Skin Friction Distribution (Upper Surface), M=0.3, a=4% 

Re=3xl06 
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Chordwise Location, x/c 
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Figure 19. Skin Friction Distribution (Lower Surface), M=0.3, a=4% 

Re=3xl06 
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Chordwise Location, x/c 
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Figure 20. 



Mach Contour, M=0.3, a=4°, Re=3xl0^ 
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Figure 21. Density Contour, M=0.3, a=4°, Re=3xl0^ 
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Convergence History, M=0.3, a=l 1.0, Re=3x 10**6 




l^npiSQ'^ XBJAJ 



Figure 22. Convergence History, M=0.3, a=ll°, Re=3xl06 
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Iteration 



M=0.3, a=l 1.0 deg, Re=3x 10**6 
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Figure 23. Pressure Distribution, M=0.3, a=ll°, Re= 3 xl 06 
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Chordwise Location, x/c 



M=0.3, a=11.0 deg, Re=3xl0**6 
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Figure 24. Skin Friction Distribution (Upper Surface), M=0.3, 

a=ll°, Re=3xl06 
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Chordwise Location, x/c 



M=0.3, a=l 1.0 deg, Re=3x 10**6 
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Figure 25. Skin Friction Distribution (Lower Surface), M=0.3, 

a=ir, Re=3xl06 
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Chordwise Location, x/c 
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Figure 26. Mach Contour, M=0.3, a=H% Re-3xl06 
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Figure 27. Density Contour, M=0.3, a=ll°, Re=3xlO^ 
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Convergence History, M=0.3, a=14.0, Re=3xl0**6 




Figure 28. Convergence History, M=0.3, a=14°, Re=3xl06 
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Iteration 



Steady State: M=0.3, Re=3x 10**6 
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Figure 29. Pressure Distribution, M=0.3, a=14“, Re=3xl06 
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Chordwise Location, x/c 



Steady State: M=0.3, Re=3x 10**6 
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Figure 30. C| vs a, M=0.3, Re= 3 xl 06 
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Steady State: M=0.3, Re=3x 10**6 
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Figure 31. Cm vs a , M=0.3, Re=3xl06 
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d. Case 4. M=0.3, Re=4xlO^ 

Figure 32 shows a comparison of the computed lift curve with the 
experimental data at a Reynolds number of four million. It can be seen that the 
computed maximum lift angle and hence the static stall angle is 13.5° degrees. 
This agrees quite well with the experimentally measured stall angle of 
McCroskey et al [Ref. 17]. The computed lift and pitching moment are in fair 
agreement for small angles of attack, but start to deviate at larger angles. 
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Figure 32. C| vs. a and Cm vs. a, M=0.3, Re=4xl06 



73 



2. Rapidly Pitching Airfoil 

Next, AGARD Case 6 from Landon [Ref. 18] was investigated using 
Ns2.f and U2diif.f. The computed results are compared to the measured data. 
The flow condition for this case corresponds to a freestream Mach number 
equal to 0.3 at a non-dimensional pitch rate equal to 0.01272 and a Reynolds 
number of 2.7 million. Flow solutions for a rapidly pitching airfoil were 
computed by pitching the airfoil about the quarter chord at a constant rate from 
zero degrees angle of attack to a final angle of attack of 15.54°. Figure 33 
shows a sketch of the motion produced by a rapidly pitching airfoil. 
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Figure 33. Rapidly Pitching Airfoil Motion 

First a steady-state solution at zero degrees angle of attack was well 
converged to 3800 iterations. Then the iteration counter was reset to zero and 
the airfoil was rapidly pitched up to the final angle. The reduced frequency k is 
given by k=a c/ 2 Uoo, where a is the pitch rate. In terms of non-dimensional 
quantities a(t)=2Mook and a(t)= ao+ai(t/ T), where ao (Aq), is the initial angle 
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of attack a=0°, ai (Ai), is the final angle of attack a=15.54°, and T is the time 
required for the airfoil to complete the pitching motion. 

a. Case 5. M=0.3, k=0.01272, Re=2.7xI06 

Figure 34 shows the comparison of the computed pressure 
distribution for the ramp-type unsteady motion at a=4.84°, for an inviscid Euler 
calculation. This Euler solution uses a 121x35 C-grid for the computation. 
Similar results were produced for ensuing angles of attack. Figures 35 through 
39 show pressure distributions at a=4.84°, a=6.72°, a=10.80°, a=12.83°, and 
a=15.54°. Surprisingly, the inviscid code predicted the suction peak more 
accurately than the N-S code for all angles of attack. Several explanations are 
offered to account for this unexpected result. First the 161x65 C-grid may not 
be fine enough to predict the magnitude of the suction peak and the N-S code 
may have too much numerical viscosity, thus in effect performing calculation at 
a lower Reynolds number. The flow at the leading edge may in fact be laminar, 
while the computed solution assumed turbulent flow everywhere. Figure 40 
shows the Mach contour and Figure 41 shows the density contour at a=10.8°. 
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Cp plot for 0012 ramp M = 0.3 a = 4.87 





Figure 34. Pressure Distribution (Inviscid), Ramp Motion, M=0.3, 
a=4.87% k=0.01272, Re=2.7xl06 
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Ramp Motion; M=0.3, K=0.01272, Re=2.7x 10**6 
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Figure 35. Pressure Distribution, Ramp Motion, M=0.3, a=4.87% 

k=0.01272, Re=2.7xl06 
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Chordwise Location, x/c 



Ramp Motion: M=0.3, K=0.01272, Re=2.7x 10**6 
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Figure 36. Pressure Distribution, Ramp Motion, M=0.3, a=6.72°, 

k=0.01272, Re=2.7xl06 
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Chordwise Location, x/c 



Ramp Motion; M=0.3, K=0.01272, Re=2.7x 10**6 
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Figure 37, Pressure Distribution, Ramp Motion, M=0,3, a=10,80°, 

k=0.01272, Re=2,7xl06 
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Chordwise Location, x/c 



Ramp Motion: M=0.3, K=0. 01272, Re=2.7xl0**6 
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Figure 38. Pressure Distribution, Ramp Motion, M=0.3, a=12.83°, 

k=0.01272, Re=2.7xl06 
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Chordwise Location, x/c 



Ramp Motion: M=0.3, K=0.01272, Re=2.7x 10**6 
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Figure 39. Pressure Distribution, Ramp Motion, M=0.3, a- 15.54 , 

k=0.01272, Re=2.7xl06 
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Figure 40. Mach Contour, Ramp Motion, M=0.3, a=10.8% 

k=0.01272, Re=2.7xl06 
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Figure 41. Density Contour, Ramp Motion, M=0.3, cc-10.8 , 

k=0.01272, Re=2.7xl06 



83 



A detail of the flow field generated by Ns2.f at a= 12.83° is shown 
in Figure 42. The flow appears to be smooth and attached (the dark line shows 
the boundary layer edge). However, magnification, as shown in Figure 43, 
reveals the start of a slight reverse flow region at approximately 10% chord. At 
a=15.54°. Figure 44, the boundary layer is much thicker with reverse flow 
profiles from approximately 10% chord to the trailing edge. 

Figures 45 and 46 present the computed Cl vs a and Cm vs a 
curves, respectively. The Navier-Stokes code underpredicts lift and pitching 
moment coefficients whereas the inviscid panel method gives significantly 
better lift predictions. However, it fails to predict the pitching moment 
coefficient for angles greater than 10°. A more comprehensive investigation of 
the dynamics of a rapidly pitching airfoil can be found in Grohsmeyer [Ref. 19]. 
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Figure 42. Velocity Field, Ramp Motion, M=0.3, a=12.83°, 

k=0.01272, Re=2.7xl06 
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Figure 43. Velocity Field (Magnified), Ramp Motion, M=0.3, 
a=12.83°, k=0.01272, Re=2.7xl06 
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Ramp Motion: M=0.3, K=0.01272, Re=2.7xl0**6 
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Figure 45. Ci vs (X, Ramp Motion, M=0.3, k=0. 01272, Re=2.7xl0^ 
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Ramp Motion: M=0.3, K=0.01272, Re=2.7xl0**6 
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Figure 46. Cm vs a , Ramp Motion, M=0.3, k=0.01272, Re=2.7xl0^ 



89 



Alpha, De 



3. Harmonic Airfoil Oscillation Using the Baldvvin-Lomax 
Turbulence Model 

Next, the code was applied to sinusoidal airfoil oscillations about the 
quarter chord point. These time periodic solutions were obtained for the second 
cycle because the results obtained for the third cycle were indistinguishable 
from ones of the second cycle. Figure 47 shows the computed angle of attack 
history. 




Figure 47. Harmonically Oscillating Airfoil 

Here, Aq is the mean angle of attack, and Ai is the amplitude. For 
each oscillatory case a steady-state solution was obtained for the minimum 
angle of attack during the cycle. The iteration counter was set to zero as the 
sinusoidal motion was time shifted half a cycle so that the start of the motion 
begins at the minimum angle of attack. This ensures a smooth transition for the 
converged steady state solution to the start of the oscillatory motion. Ten 
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thousand (10,000) iterations were computed per cycle, with a timestep equal to 
approximately 0.05 to 0.01 . 

a. Case 6. a(t)=Ao°±A j °sin(cot) 

All the computations presented in this section are computed with 
Ns2.f using the Baldwin-Lomax turbulence model. Figures 48 and 49 show the 
effect of Mach number on the lift and moment hysteresis loops for small 
amplitude oscillations about a mean angle of 12 degrees with an amplitude of 
3.0 degrees, at a Reynolds number of four million and at a reduced frequency 
of 0.1. Note that the hysteresis loops for Moo=0.4 oscillate wildly and third 
cycle computations do not coincide with results from previous cycles. For the 
Moo=0.3 computation the moment hysteresis loop indicates a stable condition. 
The effect of amplitude is shown in Figures 50 and 51. Figures 52 and 53 also 
show the effect of amplitude for a larger mean angle of attack, Ao=14.0°, at the 
same Mach number, reduced frequency, and Reynolds number, as before. Both 
solutions oscillate in a manner that is not repeatable if the motion is allowed to 
continue for additional cycles. For this comparison the maximum angle of 
attack attained during a cycle exceeded 16.5°. A comparison of the first and 
second cycle lift and moment coefficient loops is shown in Figures 54 and 55. 
Note that the second cycle lift curve starts at a slightly greater lift coefficient 
than in the first cycle. Aside from this anomaly, the loops are in good 
agreement but exhibit no instability. 

Figure 56 shows the computed lift and moment loops for small 
amplitude oscillations about a mean angle of 13 degrees with an amplitude of 
2.5 degrees, at a Reynolds number of two million, for three values of reduced 
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frequency. It can be seen that all the moment loops computed in this section 
using the Baldwin-Lomax turbulence model indicate torsional stability. 
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Figure 48. Effect of Mach Number on Unsteady Motion, 
a(t)=12°±3°sin(tot), k=0.10, Re=4xl06 
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Figure 49. Effect of Mach Number on Unsteady Motion, 
a(t)=12°±3°sin((0t), k=0.10, Re=4xl06 
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Figure 50. Effect of Amplitude on Unsteady Motion, M=0.3, 
a(t)=12°±Ai°sin(cot), k=0.10, Re=4xl06 
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Figure 51. Effect of Amplitude on Unsteady Motion, M=0.3, 
a(t)=12°±Ai°sin(0)t), k=0.10, Re=4xl0^ 
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Figure 52. Effect of Amplitude on Unsteady Motion, M=0.3, 
a(t)=14°±Ai°sin(cot), k=0.10, Re=4xl06 
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Figure 53. Effect of Amplitude on Unsteady Motion, M=0.3, 
a(t)=14°±Ai°sin((0t), k=0.10, Re=4xl06 
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Figure 54. Effect of 3nd 2*^^ Cycle on Unsteady Motion, M=0.3, 
a(t)=13“± 2°sin(0)t), k=0.10, Re=4xl06 
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Figure 55. Effect of l^t and 2*^^ Cycle on Unsteady Motion, M=0.3, 
a(t)=13±2°sin(CDt), k=0.10, Re=4xl06 
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Figure 56. Effect of Reduced Frequency on Unsteady Motion 
M=0.3, a(t)=13±2.5°sin((Ot), Re=2xl06 
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b. Case 7. a(t)=13'‘±2.5’’sin(a)t), k=0.2, Re=4xlO^ 

In an effort to determine if the flow field was in fact predicting a 
trailing edge vortex that would shed during the oscillation cycle as seen in the 
experimental results, instantaneous particle traces and the corresponding 
velocity fields were plotted for the downstroke at a Mach number of 0.3, a 
mean angle of attack of 13 degrees, an amplitude of 2.5 degrees, a reduced 
frequency of 0.1 and a Reynolds number of four million. Figures 57 through 61 
show the downstroke portion of the cycle beginning at a=15.3° and ending at 
a=14.5°. A recirculatory region is seen in Figure 57 which grows slightly, as 
seen in Figure 58. This region then appears to diminish in size and intensity. 
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Figure 57. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.3° (Downstroke), 
a(t) = 13°±2.5‘’sin(cot), k=0.20, Re=4xl06 
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Figure 58. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.1° (Downstroke), 
a(t) = 13°±2.5°sin(cot), k=0.20, Re=4xl0^ 
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Figure 59. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.9° (Downstroke), 
a(t)=13°±2.5°sin(cot), k=0.20, Re=4xl06 
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Figure 60. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.7° (Downstroke), 
a(t) = 13°±2.5“sin(0)t), k=0.20, Re=4xl06 
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Figure 61. Instantaneous Particle Trace and Velocity Field, 



Oscillatory Motion, M=0.3, a=14.5° (Downstroke), 
a(t)=13°±2.5°sin(cot), k=0.20, Re=4xl06 
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c. Cases. a(t)=13’'±2.5’‘sin((Ot), k=0.1, Re=2xl0^ 

The plots generated for Case 7 did show a recirculatory region, 
but a full understanding of its development and structure could not be gained 
from only a portion of the oscillatory cycle. Therefore in Figures 62 through 80 
a complete cycle is shown for the same flow condition as in Case 7 with only 
the reduced frequency being increased from 0.1 to 0.2 and Reynolds Number 
decreased from four million to two million. A recirculatory region is seen to 
form at a=14.9° on the upstroke, which continues to grow and intensify. This 
recirculatory region reaches its largest value at a=15.3° on the downstroke and 
then diminishes and finally the flow becomes completely reattached. It is 
important to note that even for this relatively fast oscillation the recirculatory 
region did not shed into the wake. 
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Figure 62. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=12° (Upstroke), a(t)=13°±2.5°sin(o)t), 

k=0.10, Re=2xl06 
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Figure 63. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=13° (Upstroke), a(t)=13°±2.5°sin(o)t), 

k=0.10, Re=2xl06 



110 





normal I -20. pressure 
=J3.D, fl=2.5. k=O.I. B-L t»r Miusiroke) 



^C^TOUR LEUFLS 
n sRiifin ^ 



0.300 MACH 
SfEG ALPHA 
2 . OTJh-I^ *GR e 
1 . 3-4x 10*^>tj_ME 

grTD^- 




Figure 64. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=13.5° (Upstroke), a(t)=13°±2.5°sin(cot), 

k=0.10, Re=2xl06 
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Figure 65. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14° (Upstroke), a(t)=13°±2.5°sin(o)t), 

k=0.10, Re=2xl06 
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Figure 66. Velocity Field, Oscillatory Motion, M=0.3, a=14.5' 
(Upstroke), a(t)=13°±2.5°sin(cot), k=0.10, Re= 2 xl 06 
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Figure 67. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.7° (Upstroke), a(t)=13°±2.5°sin(o)t), 

k=0.10, Re=2xl06 
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Figure 68. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.9° (Upstroke), a(t)=13°±2.5°sin((Ot), 

k=0.10, Re=2xl06 
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Figure 69. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.1° (Upstroke), a(t)=13°±2.5''sin((Ot), 

k=0.10, Re=2xl06 
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Figure 70. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.3° (Upstroke), a(t)=13°±2.5'sin((0t), 

k=0.10, Re=2xl06 
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Figure 71. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.5' (Upstroke), a(t)=13'’±2.5°sin((0t), 

k=0.10, Re=2xl06 
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Figure 72. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.3° (Downstroke), 
a(t)=13°±2.5°sin((ot), k=0.10, Re=2xl0^ 
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Figure 73. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=15.1° (Downstroke), 
a(t)=13“±2.5°sin(0)t), k=0.10, Re=2xl06 
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Figure 74. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.9° (Downstroke), 
a(t)=13'±2.5°sin(cot), k=0.10, Re=2xl06 
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Figure 75. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.7° (Downstroke), 
a(t)=13°±2.5°sin((0t), k=0.10, Re= 2 xl 06 
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Figure 76. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14.5° (Downstroke), 
a(t)=13°±2.5°sin((0t), k=0.10, Re=2xl06 
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Figure 77. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=14° (Downstroke), 
a(t)=13“±2.5°sin((0t), k=0.10, Re=2xl06 
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Figure 78. Velocity Field, Oscillatory Motion, M=0.3, a=13.5° 
(Downstroke), a(t)=13°±2.5°sin(©t), k=0.10, Re=2xl0^ 
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Figure 79. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=13° (Downstroke), 
a(t)=13°±2.5°sin(0)t), k=0.10, Re=2xl06 
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Figure 80. Instantaneous Particle Trace and Velocity Field, 
Oscillatory Motion, M=0.3, a=12° (Downstroke), 
a(t)=13°±2.5°sin(cot), k=0.10, Re=2xl06 
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4. Effect of Turbulence Modeling 

a. Case 9. a(t)=9‘‘±5°sin((ot), k=0.2, Re=4xl0^ 

Figures 81 and 82 show a comparison of the lift and moment 
hysteresis loops using the Baldwin-Lomax turbulence model for an oscillation of 
5 degrees amplitude about a mean angle of attack of 9.0 degrees with 
McCroskey’s experimental data for a Mach number of 0.3, a Reynolds number 
of four million, and a reduced frequency k=0.2. It is seen that there is a 
substantial difference between the computed and the experimental lift 
hysteresis loops. This difference becomes even more pronounced when the 
computed and experimental moment hysteresis loops are compared with each 
other. In an effort to understand the failure of Ns2.f to predict the experimental 
pitching moment hysteresis loop, the pressure distributions for several upstroke 
and downstroke angles of attack were compared to measured data from 
Reference 1. Figures 83 through 87 show pressure distributions for a=11.9° on 
the upstroke portion of the cycle; and a=13.7°, a=13.0°, a=10.5°, and a=8.9° on 
the down stroke cycle. U2diif.f overpredicts the suction peak at a=11.9° on the 
upstroke and Ns2.f underpredicts the peak. Even though the values are off, the 
general shape of the distribution is well predicted. However, on the downstroke 
at a=13.7°, a significant difference between the computed and the measured 
distributions is observed. The measured pressure distribution shows a lower 
suction near the leading edge but higher suction near the trailing edge. Neither 
of these changes are predicted by Ns2.f. As the oscillation continues on the 
downstroke the measured distribution tends to flatten out creating a “plateau”. 
Just after a=8.9° on the downstroke, the flow begins to reattach and the positive 
pressure difference at the trailing edge becomes smaller. The failure of Ns2.f to 
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predict the experimental loop results from the inability of Ns2.f to accurately 
compute the pressure distributions. Since the pitching moment is an integrated 
quantity, the size of the area under the pressure distribution aft of the quarter 
chord must be greater than the area forward of the quarter chord to compute a 
negative or downward pitching moment. Ns2.f clearly does not predict the 
experimental pressure distributions during the downstroke. 
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Oscil. Motion: M=0.3, Ao=9, Al=5, k=.2, Re=4xl0**6 




juaioyjaoo yn 



Figure 81. C| vs a , Oscillatory Motion, M=0.3, a(t)=9°±5°sin(cot), 

k=0.2, Re=4xl06 
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Oscil. Motion: M=0.3, Ao=9, Al=5, k=.2, Re=4xl0**6 




I I 



3U9pyj903 5U91LIOJAI 

Figure 82. Cm vs a , Oscillatory Motion, M=0.3, a(t)=9°±5°sin(0)t), 

k=0.2, Re=4xl0^ 
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Oscil. Motion: M=0.3, Ao=9, A 1=5, k=.2, Re=4x 10**6 
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Figure 83. Pressure Distribution, Oscillatory Motion, M=0.3, 
0=11.9“ (Upstroke), a(t)=9“±5“sin(o)t), k=0.2 Re= 4 xl 06 
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Oscil. Motion: M=0.3, Ao=9, Al=5, k=.2, Re=4xl0**6 




p 



oo 

o 



kO 

o 



o 



o 



o 

o 



5U9ioyj903 9jnss9Jj - 



Figure 84. Pressure Distribution, Oscillatory Motion, M=0.3, 
a=13.7° (Downstroke), a(t)=9°±5°sin(o)t), k=0.2, Re=4xl0^ 
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Oscil. Motion: M=0.3, Ao=9, A 1=5, k=.2, Re=4x 10**6 
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Figure 85. Pressure Distribution, Oscillatory Motion, M=0.3, 
a=13.0° (Downstroke), a{t)=9°±5°sin(o)t), k=0.2, Re=4xl0^ 
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Oscil. Motion: M=0.3, Ao=9, A 1=5, k=.2, Re=4x 10**6 
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Figure 86. Pressure Distribution, Oscillatory Motion, M=0.3, 
0=10.5” (Downstroke), a(t)=9°±5°sin((ot), k=0.2, Re= 4 xl 06 
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Chordwise Location, x/c 



Oscil. Motion: M=0.3, Ao=9, Al=5, k=.2, Re=4xl0**6 
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Figure 87. Pressure Distribution, Oscillatory Motion, M-0.3, a-8.9 
(Downstroke), a(t)=9°±5°sin(©t), k=0.2, Re=4xl0^ 
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b. Case 10. a(t)=10°±5.5^sin((Ot), k=0.1, Re=4xlO^ 

The measured loop consists of two subloops, a clockwise 
“unstable” subloop and a counterclockwise “stable” subloop. It is seen that the 
Navier-Stokes calculations in combination with the Baldwin-Lomax turbulence 
model fails to predict the destabilizing clockwise pitching moment loop. For this 
reason the sensitivity of the computed loops to different turbulence models was 
studied next. Figures 88 and 89 show the computed lift and pitching moment 
loops for oscillation about a slightly higher mean angle of 10 degrees at an 
amplitude of 5.5 degrees. The Mach and Reynolds numbers are again 0.3 and 
four million, and k was decreased to 0.1 to minimize the effects of reduced 
frequency. It can be seen that the Baldwin-Lomax, Johnson-King, and the 
RNG turbulence models produce significantly different hysteresis loops. The 
RNG model produces relatively good agreement with the measured lift 
hysteresis loop but fails again to predict the destabilizing moment subloop. 
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Figure 88. Ci vs a , Oscillatory Motion, M=0.3, 
a(t)=10°±5.5°sin(o)t), k=0.10, Re=4xl06 
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Figure 89. Cm vs a , Oscillatory Motion, M=0.3, 
a(t)=10°±5.5°sin((Ot), k=0.10, Re=4xl06 
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c. Case 11. a(t)=ll'‘±5'‘sin((0t), k-0.1, Re=4xlO^ 

A solution was also computed for oscillations about a mean angle 
of 1 1 degrees at an amplitude of 5 degrees, the other parameters left 
unchanged. Figures 90 and 91 show the computed lift and moment hysteresis 
loops. As before, the calculations fail to predict the expected destabilizing 
moment subloop. The reason for this failure can be better appreciated by 
visualizing the separated flow structures computed with these turbulence 
models. Figures 92 and 93 show the computed flowfield with particle traces 
and the velocity field. Figure 92 corresponds to the instant when the airfoil 
oscillates through an incidence of 15 degrees during the downstroke. Different 
turbulence models produce substantially different recirculatory flow patterns on 
the upper surface near the trailing edge. The Baldwin-Lomax model produced 
the smallest recirculatory region while the RNG model produced the largest. In 
addition. Figure 93 (a=15°) shows the relative magnitudes of the velocity 
vectors. Note, near the surface at the trailing edge, the magnitude of the 
velocity vectors in the reverse flow region calculated by the Baldwin-Lomax 
model are much smaller than the region calculated by the Johnson-King model 
or RNG model. As the cycle continues. Figures 94 and 95 (a=14°) show the 
recirculation region to shrink slightly and loose intensity for the Baldwin-Lomax 
model, while the regions as calculated by the Johnson-King and RNG continue 
to expand and intensify. As noted in Case 8 the recirculatory flow region 
grows and then vanishes again during the course of the oscillation. Even 
though the structure at the trailing edge resembles a vortex it does not shed 
from the trailing edge. The occurrence of the destabilizing moment loop 
appears to be closely associated with the vortex shedding from the trailing 
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edge, an event which occurs soon after the static angle is substantially 
exceeded during part of the cycle. Although the static stall angle is significantly 
exceeded in the present calculations the numerical procedure along with the 
turbulence modeling used in the calculations appears to be unable to produce 
the anticipated shedding process. 
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Figure 90. C| vs a , Oscillatory Motion, M=0.3, a(t)=ll°±5°sin(©t), 

k=0.10, Re=4xl06 
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Figure 91. Cm vs a , Oscillatory Motion, M=0.3, a(t)=ir±5°sin(o)t), 

k=0.10, Re=4xl06 
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Angle of attack, deg 



EFFECT OF TURBULENCE MODELING ON THE COMPUTED 
TRAILING EDGE FLOWFIELD 

= 0.3, a (t) = H'’ + 5- sin (wt), k = 0.1, Re = 4 x 10® 
a =15 DOWNSTROKE 




(a) Baldwin-Lomax 




(b) Johnson-King 




(c) RNG 



Figure 92. Instantaneous Particle Trace, Oscillatory Motion, M=0.3, 
a=15° (Downstroke), a(t)=ll°±5°sin(cot), k=0.10, Re=4xl0^ 
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EFFECT OF TURBULENCE MODELING ON THE COMPUTED 
TRAILING EDGE FLOWFIELD 

Moo = 0.3, a (t) = 11° + 5° sin (cot), k = 0.1, Re = 4 x 10^ 
a = 15’ DOWNSTROKE 




(a) Baldwin-Lomax 





Figure 93. Velocity Field, Oscillatory Motion, M=0.3, a=15° 
(Downstroke), a(t)=ll°±5°sin(0)t), k=0.10, Re=4xl06 
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EFFECT OF TURBULENCE MODELING ON THE COMPUTED 
TRAILING EDGE FLOWFIELD 

M = 0.3, a (!) = ir + 5° sin (wt), k = 0.1, R© = 4 x 106 
a = 14® DOWNSTROKE 




(a) Baldwin-Lomax 




(b) Johnson-King 




(c) RNG 



Figure 94. Instantaneous Particle Trace, Oscillatory Motion, M=0.3, 
a=14° (Downstroke), a(t)=ll°±5°sin(o)t), k=0.10, Re=4xl0^ 
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EFFECT OF TURBULENCE MODELING ON THE COMPUTED 
TRAILING EDGE FLOWFIELD 

M = 0.3, a (t) = ir + 5’ sin (cat), k = 0.1, Re = 4 x 10^ 
a = 14” DOWNSTROKE 






Figure 95. Velocity Field, Oscillatory 
(Downstroke), a(t)=ll°±5°sin(cot), 



Motion, 

k=0.10. 



M=0.3, a=14° 
Re=4xl0^ 
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V. CONCLUSIONS AND RECOMMENDATIONS 



In the preceding investigation a computationally efficient, fully factorized, 
two-dimensional Navier-Stokes flow solver was utilized to predict steady and 
unsteady flow solutions about a NACA 0012 airfoil. Comparisons between a 
steady viscous/inviscid interaction method code using the Cebeci-Smith 
turbulence model and the Navier-Stokes code using the Baldwin-Lomax 
turbulence model show close agreement up to just prior to the static stall angle 
of attack. After this point, the viscous/inviscid interaction method code fails to 
accurately model the flow separation. 

When applied to the unsteady problem of airfoil stall flutter in compressible 
flow, the Navier-Stokes code shows that the modeling of the recirculatory flow 
region with current turbulence models fails to capture the essential physics 
which governs the onset of stall flutter. Comparisons of the computed results 
with available experimental data indicates that even though the lift response is 
fairly well predicted, the computation of the pitching moment hysteresis loops is 
very sensitive to turbulence modeling. Results computed with several current 
models are in good agreement whenever the steady stall angle is exceeded only 
slightly. However, they fail to capture a vortex shedding process that may 
contribute to the onset of stall flutter. 

Therefore, further detailed studies of improved numerical schemes and 
turbulence models as well as viscous/inviscid interaction approaches are 
required to improve the prediction of the unsteady flow separation and vortex 
shedding phenomenon. Further computational work with the full Navier-Stokes 
equations instead of applying the thin-layer approximation is recommended. In 
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addition, local grid refinement studies that focus on the the leading and trailing 
edge upper surface may lead to the accurate prediction of the extremely 
complex vortex development and shedding process. 
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APPENDIX A - FLUX JACOBIAN MATRIX 



The Jacobian matrices = dE ! dq and B = <3F/(^ are given by A or B = 
kQl + kiA + k2B, where A = dE/3q and B = 3F/8q are the usual Jacobian 
matrices of the Cartesian flux vectors. The Aor B matrix is 



r ko 



k, 



kj 0 1 




-u(kjU +kjv) 
+ 

-v(k;U +kjV) 

+ k^^^ 



-(7-2)k,u 

+ k()+ kjU+ kjV 

k.v 

-fy-OkjU 



-(y-l)k,v 

+ kjU 



(y-i)k, 



-(r-2)k2v+ 
ko +k,u +kjV 



|(k,u+k2v)* |y(e/p)-(/>2|^,^ _ [^(e/p)_p2j,^^_ y(k,u + k2v)| 

[l“7(s/p)+2^] (y - i)(k,u + k2v)u (y- l)(k,u + k2v)v + k^ J 



where (j)2 = 0.5(y - l)(u2 + v2), kQ=^t, ki=^x, k2=^z for A, and ko=Ct, ki=Cx, 
k2=Cz for B. 
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APPENDIX B - PITCHING MOMENT DERIVATION 



Figure 96 shows the nomenclature for the integration of the pressure and 
shear stress distributions over a two dimensional body surface as generated by 
the Navier-Stokes code Ns2.f. For unsteady motion the airfoil is rotated with 
respect to a fixed inertial frame of reference (x,z; Laboratory frame of 
Reference). Note, for unsteady motion, as the airfoil rotates through an angle of 
attack, all grid (i=l...Imax» k=l...Kmax) iire rotated about the desired pivot point 
with respect to the fixed x,z coordinate system. As the solution is advanced in 
time the flow quantities for this new grid location are updated for each point. 



I e 



/ 





k 





T.E. upper, i=131 



T.E. lower, i=3 1 



S 



Ax=x(i+1) - x(i) 



X - coordinate 



Figure 96. Pitching Moment Derivation 
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PITCHING MOMENT DERIVATION GENERATED BY MATHCAD 



1 . Read tabulated data generated from Ns2.f. Data computed at a=0° and a Re=4x 10**6. 

Data : = READPRN(up) Data := READPRN(low) i := 1 ..50 

u 1 



u 



u 





<0> 




<0> 


X 


:= Data 


X 


:= Data 


u 


u 


1 


1 




<1> 




<1> 


z 


:= Data 


z 


;= Data 


u 


u 


1 


1 




<2> 




<2> 


cp 


:= Data 


cp 


:= Data 


u 


u 


1 


1 




<3> 




<3> 


T 


:= Data 


T 


:= Data 


u 


u 


1 


1 


and Az for upper and lower surface. 






i X 


- X 


A X . 1 


: = X . 1 - X . 1 


u 


u 








i i-1 


i 


i ( 


: Z 


- z 


A z . 1 


:= z.l - z.l 


U 


u 








i i-1 


i 


i ( 



3. Compute segment length As for upper and lower surface. 



6 s 



u 























_ 






2 




2 








2 




2 


= 




Ax 


+ 


A z 


As. 1 := 






A X . 1 


+ 


’a z.l 








u 




u 
















> 










i > 















4. Compute angle, 6, relative to perpendicular for upper and lower surface. 





A z 








A z 




u 








1 


IT 




i 


o 

II 

CD 


1T 


+ atari 


i 




+ atari 








180 




Ax 


1 


.180. 




A X 




u 


i 






1 












i- 
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Computed the x and z moment arms for the upper and lower surface. 



X + X 

u u 

i i-1 



z + z 
u u 

i i-1 



XX 



z z 



u 



X + X 

1 1 

i i-1 



z + z 
1 1 

i i-1 



XX 



z z 



. Computepitching moment coefficient about quarter chord by applying equation 1.11 from [Ref. 
, p. 17J. 



Cm 



25 



i 



f-cp "COS 


‘a 


+ T ‘sin 


‘a 




u 


u 


u 


L i 


i_ 


i 


i_ 



• XX 



u 



‘ A s 



fcp ’sin 


"a 


+ T ’COS 


’a 




’ zz 




u 


u 


u 




u 


L i 


i_ 


i 


i_ 




i_ 


Cp ’COS 


*a 


+ T ‘ sin 


'a 


- 


‘ XX 


1 


1 


1 


1 




1 


i 


i_ 


i 


i_ 




i 



u 



• A s 



-cp “Sin 


‘a 


+ T “COS 


"a 


“ z z 


1 


1 


1 


1 


1 


i 


i^ 


i 


i. 


i_ 



Cm = 0.005 
25 

. The computed pitching moment coefficient from Ns2d.f (Computed on a Cray Super-Computer) 
/as Cm=0.(X)06, as compared to the calculations made on a Macintosh SE/30 where Cm=0.()05. 
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Upper Surface Data: 



X 


z 


-Cp 


Tau 


-2.5e-l 


-1.6e-3 


-1.0319e+0 


3e-3 


-2.497e-l 


3e-3 -9. 


841e-l -7.4e- 


•3 


-2.479e-l 


7.6e-3 


-6.798e-l 


-1.5e-2 


-2.447e-l 


1.2e-2 


-4.091e-l 


-1.64e-2 


-2.404e-l 


1. 65e-2 


-1.766e-l 


-1.54e-2 


-2.346e-l 


2.09e-2 


9.93e-2 


-1.19e-2 


-2.273e-l 


2.49e-2 


2.286e-l 


-8.4e-3 


-2.187e-l 


2.88e-2 


2.769e-l 


-6.3e-3 


-2.088e-l 


3.26e-2 


3.316e-l 


-4 . 9e-3 


-1.975e-l 


3.61e-2 


3.715e-l 


-3.8e-3 


-1.85e-l 


3.94e-2 


3.951e-l 


-3e-3 


-1.713e-l 


4.26e-2 


4,105e-l 


-2.4e-3 


-1.564e-l 


4.55e-2 


4.199e-l 


-1.9e-3 


-1.403e-l 


4.81e-2 


4.245e-l 


-1.5e-3 


-1.23e-l 


5.05e-2 


4,248e-l 


-1.2e-3 


-1.047e-l 


5.26e-2 


4 .221e-l 


-le-3 


-8.54e-2 


5.44e-2 


4 .168e-l 


-7e-4 


-6.51e-2 


5.6e-2 


4.094e-l 


-6e-4 


-4.38e-2 


5.73e-2 


3.996e-l 


-4e-4 


-2.17e-2 


5.82e-2 


3.904e-l 


-3e-4 


1.3e-3 


5.89e-2 


3.778e-l 


-2e-4 


2.5e-2 


5.93e-2 


3.628e-l 


-le-4 


4.95e-2 


5.94e-2 


3.494e-l 


Oe+0 


7.45e-2 


5.92e-2 


3.343e-l 


le-4 


l,002e-l 


5.88e-2 


3.184e-l 


le-4 


1.264e-l 


5.81e-2 


3.031e-l 


2e-4 


1.53e-l 


5.72e-2 


2.875e-l 


2e-4 


1.8e-l 


5.6e-2 


2.71e-l 


2e-4 


2.073e-l 


5.46e-2 


2.544e-l 


3e-4 


2.348e-l 


5.3e-2 


2.378e-l 


3e-4 


2.626e-l 


5.13e-2 


2.211e-l 


3e-4 


2.904e-l 


4.93e-2 


2.045e-l 


3e-4 


3.182e-l 


4.72e-2 


1.884e-l 


3e-4 


3.461e-l 


4.5e-2 


1.716e-l 


4e-4 


3.738e-l 


4.27e-2 


1.559e-l 


4e-4 


4.014e-l 


4.02e-2 


l,421e-l 


4e-4 


4.287e-l 


3.77e-2 


l,252e-l 


4e-4 


4.557e-l 


3.5e-2 


1.067e-l 


4e-4 


4.823e-l 


3.23e-2 


9.05e-2 


4e-4 


5.085e-l 


2.96e-2 


7.46e-2 


4e-4 


5.342e-l 


2.68e-2 


5. 68e-2 


4e-4 


5.593e-l 


2.4e-2 


3.86e-2 


4e-4 


5.838e-l 


2.12e-2 


2.06e-2 


4e-4 


6.076e-l 


1. 84e-2 


1. 3e-3 


3e-4 


6.306e-l 


1.56e-2 


-2.05e-2 


3e-4 


6.528e-l 


1.29e-2 


-4.42e-2 


3e-4 


6.742e-l 


1. Ole-2 


-7.02e-2 


3e-4 


6.946e-l 


7.4e-3 


-9.2e-2 


2e-4 


7.141e-l 


4.7e-3 


-1.441e-l 


le-4 


7,325e-l 


2.2e-3 


-1.088e-l 


Oe+0 


7.499e-l 


Oe+0 -2. 


34e-l Oe+0 
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Lower Surface Data: 



X 


z 


-Co 


Tau 


-2.5e-l 


-1. 6e-3 


-1.0319e+0 


3e-3 


-2.486e-l 


-6.2e-3 


-7.832e-l 


4 . 9e-3 


-2.458e-l 


-1.06e-2 


-4.815e-l 


2.3e-3 


-2.419e-l 


-1.51e-2 


-2.554e-l 


-2e-4 


-2.366e-l 


-1.96e-2 


2.4e-2 


-2.2e-3 


-2,298e-l 


-2.37e-2 


2.063e-l 


-2.8e-3 


-2.215e-l 


-2.76e-2 


2.624e-l 


-2.7e-3 


-2.12e-l 


-3.14e-2 


3.142e-l 


-2.5e-3 


-2,012e-l 


-3.5e-2 


3.619e-l 


-2.2e-3 


-1.891e-l 


-3.84e-2 


3.892e-l 


-1.9e-3 


-1.758e-l 


-4.16e-2 


4.067e-l 


-1.7e-3 


-1.613e-l 


-4.45e-2 


4.178e-l 


-1.4e-3 


-1.456e-l 


-4.73e-2 


4.239e-l 


-1.2e-3 


-1.288e-l 


-4.97e-2 


4.256e-l 


-le-3 


-1.109e-l 


-5.19e-2 


4.237e-l 


-8e-4 


-9.19e-2 


-5.39e-2 


4.l91e-l 


-7e-4 


-7.2e-2 


-5.55e-2 


4.126e-l 


-5e-4 


-5.11e-2 


-5.69e-2 


4.035e-l 


-4e-4 


-2.94e-2 


-5.79e-2 


3.932e-l 


-3e-4 


-6.8e-3 


-5.87e-2 


3.84e-l 


-2e-4 


1 . 66e-2 


-5.92e-2 


3.681e-l 


-le-4 


4.06e-2 


-5. 94e-2 


3.527e-l 


Oe+0 


6.53e-2 


-5. 93e-2 


3.419e-l 


le-4 


9.06e-2 


-5.9e-2 


3.249e-l 


le-4 


l.l64e-l 


-5.84e-2 


3.084e-l 


2e-4 


1.427e-l 


-5.76e-2 


2.941e-l 


2e-4 


1.694e-l 


-5.65e-2 


2.777e-l 


3e-4 


1.963e-l 


-5.52e-2 


2.612e-l 


3e-4 


2.236e-l 


-5.37e-2 


2.448e-l 


3e-4 


2.51e-l 


-5.2e-2 


2.283e-l 


4e-4 


2.785e-l 


-5.02e-2 


2.118e-l 


4e-4 


3.061e-l 


-4.82e-2 


1.957e-l 


4e-4 


3 . 336e-l 


-4. 6e-2 


1.794e-l 


4e-4 


3.61e-l 


-4.38e-2 


1.631e-l 


1 

<D 


3.883e-l 


-4.14e-2 


1.49e-l 


4e-4 


4 .154e-l 


-3.89e-2 


1.342e-l 


4e-4 


4.421e-l 


-3.64e-2 


1.16e-l 


4e-4 


4.685e-l 


-3.37e-2 


9.88e-2 


5e-4 


4.945e-l 


-3.11e-2 


8.37e-2 


5e-4 


5.199e-l 


-2.84e-2 


6.73e-2 


5e-4 


5.448e-l 


-2.57e-2 


4. 95e-2 


5e-4 


5.69e-l 


-2.29e-2 


3.18e-2 


4e-4 


5.926e-l 


-2.02e-2 


1.43e-2 


4e-4 


6.155e-l 


-1.75e-2 


-5.3e-3 


4e-4 


6.375e-l 


-1.48e-2 


-2. 69e-2 


4e-4 


6.587e-l 


-1.21e-2 


-5.06e-2 


4e-4 


6.789e-l 


-9.5e-3 


-1 .Ale-2 


4e-4 


6.982e-l 


-6.9e-3 


-9.64e-2 


3e-4 


7.165e-l 


-4.4e-3 


-1.491e-l 


2e-4 


7.338e-l 


-2.1e-3 


-l.lOle-1 


le-4 


7.499e-l 


Oe+0 -2. 


. 34e-l Oe+0 
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APPENDIX C 



USING THE PROGRAMS 



A. DIRECTIONS FOR THE EXECUTION OF NS2.F 
1. Steady State Case 

In order to generate a Cp, Q, or Cm plot from Ns2.f there are many steps, using 
six different programs that must be taken: (1.) grape (grid generation program), (2.) rotgr 
(grid rotation program), (3.) Ns2.f (Navier-Stokes code), (4.)vi (editor), (5.) pIot3d 
(flow visualization program), and (6.)xyplot (a simple x-y plot program). 

A. ) Ns2.f is a fortran code and will take hours to run even on the Stardent mini- 
super computer. To run this code, an executable file (ns2), an input file ns2.in (actual 
name of input file is users choice) and a grid file FOROl l.DAT (Stardent) or fort. 1 1 (on 
other computers) must be in the directory. 

B. ) To obtain a grid, a grid generation program such as grape must first be 
used. Grape output places the leading edge of the grid at the reference axis at 0.0° AOA. 
Now, the grid must be rotated to the desired AOA by using the program rotgr. The input 
file (output from grape) for rotgr must be FOR001.DAT and the output is FOR002.DAT. 
The output from rotgr (FOR002.DAT) must be renamed to FOROl l.DAT for use in ns2. 

C. ) Next edit the input file to reflect the desired flow conditions. Any editor will 
do, however vi is perhaps the most common editor and it is found on most unix based 
machines. Enter Mach No. and AOA, check smoothing factors, set time step, Courant No. 
and initial no. of iterations (small value- for first run). Check restart false for first run and 
make sure No. point of your grid match input values. 

D. ) The code is executed as follows: ns2 <ns2.in> ns2.out 

The first ten iterations should take 1-10 minutes, depending on the type of 
computer used. Several output files will be generated: 

TABLE 3. NS2.F OUPUT 



OUTPUT 


COMMENTS 


ns2.out: 


convergence data 


FOR004.DAT: 


Q data 


FOR007.DAT: 


residuals 


FOR008.DAT: 


Cd, Cm data 


FOR009.DAT: 


Cp data (to be used with xyplot) 


FOR031.DAT: 


Flow data (to be used with plot3d) 



E.) Re-run the code for a greater number of iterations to achieve a 
converged solution. When the change in the max residual (FOR007.DAT) has dropped 
10'^, the solution has converged. Ensure restart is set to true. The CFL can also be 
adjusted to obtain a converged solution more rapidly. 



158 



4 



a. Sample Input to Ns2.f 



MACH, 


ALFAO, 


ALFAl, 


ALFARE, 


REDFRE, REYNOLDS 


0.300, 


0.00, 


0.0, 


0.00, 


0.00, 2.70 


ED2X, 


ED2Y, 


ED4X, 


ED4Y, 


ED 


0.00, 


0.00, 


0.030, 


0.030, 


2.0 


DT, 


COUR, 


NITER, 


NEWTIT 




0.0020, 


2100, 


50, 


2 




RSTRT, 


OSCIL, 


RAMP, 


NPER, 


TSHIFT 


true. 


false. 


false. 


10000, 


-0.5 


TIMEAC, 


IMPLBC, 


EXPLBC, 


CIRCOR 




1, 


false. 


true. 


false 




Vise, 


BLTM, 


JKTM, 


RNGTM 




true. 


true. 


false. 


false 




I TEL, 


ITEU 








31, 


131 








lAl, 


IA2, IA3 


, IA4, 


IA5, IA6, 


IA7, IA8, IA9, lAlO 


1250, 


1350,1450 


,1550,1570,1590, 


1610,1630,1630, 1805 



false 

UNSTST (If true Time = 0.,Set TRUE only when unsteady motion starts from steady 
steady state, TRUE for steady-state ok, but for unsteady restarts must be FALSE 
for proper recording of unsteady motion) 



Mach 

AlfaO 

Alfal 

Alfare 

Redfre 

Reynolds 



Free stream Mach number 

Angle of attack, also mean angle of attack for unsteady 
Amplitude of Oscillatory motion 
Reference angle 

Reduced frequency k = omega c / 2U 
Reynolds number Re = cU/n 



ED2x 

ED2z 

ED4x 

ED4z 

ED 

ISPEC 



X-direction 2nd order explicit smoothing ( 

0.25 < 

Z-direction 2nd order explicit smoothing ( 

0.10 < 

X-direction 4nd order explicit smoothing ( 

Z-direction 4nd order explicit smoothing ( 

Scaling of Implicit smoothing 
Spectral radious parameter 



e2x = 0.00 
e2x <0.50 
e2z = 0.00 
e2z < 0.20 
e4x =0.03 
e4x = 0.05 
e4z = 0.03 
e4z = 0.05 



subsonic, 

transonic) 

subsonic, 

transonic) 

subsonic, 

transonic) 

subsonic, 

transonic) 



Dt 

Cour 

Niter 

Newtit 



Time step 

Courant number Cu = dt ^ L_jnax 

Number of Iteration in this run 

Newton subiteration within each timestep 



RSTRT 

OSCIL 

RAMP 

NPER 

TSHIFT 



Restart 

Oscillatory motion A(t) = AO + A1 * sin ( k ^ M ^ t ) 

Ramp motion 

Number of time steps in one period of oscillation, dt=T/Nper 
Time shift in radians to start oscilation for any a(t) 



TIMEAC 

IMPLBC 

EXPLBC 



Time accurate Tacc=l and for Jacobian Scaled Dt, Tacc=0 
Implicit wall BC Treatment 
Explicit wall BC Treatment 
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Vise 

TURBL 

JKTM 

RNGTM 


: Only laminar viscosity 

: Baldwin-Lomax eddy viscosity 

: Johnson-King eddy viscosity 

: RNG eddy viscosity 


ITEL, 


ITEU : Lower and Upper trailing edge I locations 



lAl etc. : Write out (IAn/100) angles of attack in units 

61 - 70 

Read grid from unit fort. 11 and the flow from unit fort. 31 
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b. Sample Output from Ns2.f 

Grid dimensions : 161x64 
Mach = 0.30000 

Ao = 0.00000 A1 = 0.00000 k = 0.01272 

Re = 0.2700E+07 Dt = 0.00200 Cu=2100 . 00000 



Dimpl = 
D2x = 
D4x = 


2.00000 

0.00000 D2z = 0.00000 

0.03000 D4z = 0.03000 


Restart = 

Oscil 

Ramp 


= T 

= F 

= F 


Itel = 


31 Iteu = 131 He = 81 


Timeac 
BC impl 
BC expl 


= T 

= F 

= T 


Cour = 
L max = 
Dt 


2100. 

232096.4825513 

0.00904796 



It drmax dumax dvmax demax ir kr 

3000 0.573103E-07 0.371716E-07 0.563395E-07 0.129865E-06 140 63 

Ni= 1 0.903247E-08 0.711341E-08 0.803250E-08 0.219128E-07 

Alfa(t) = 0.000 Cl = 0.007196 Cd = 0.003260 Cm = -0.000692 
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2. Unsteady Case 

The execution of an unsteady case is essentially the same as for the 
steady-state case. 

A. First, a well converged steady-state solution must be obtained. 
This solution then becomes the starting point for the unsteady solution. 

B. The restart is left set to true, however, to reset the time counter 
to time=0, the UNSTST variable must be set to true for several iterations and 
then back to false for normal time counting. 

C. As the unsteady motion begins to march in time Flow data and 
Grid position data may be stored at user determined angles of attack. To store 
flow data at a=12.5°, a=13.5°, a=14.5°, etc... enter the data shown in Table 4. 
into the input file. 



TABLE 4. ANGLE OF ATTACK INPUT 



lAl 


1A2 


1A3 


1A4 


1A5 


1A6 


1A7 


IA8 


1A9 


lAlO 


1250 


1350 


1450 


1550 


1570 


1590 


1610 


1630 


1630 


1805 



D. The Flow data and the Grid data at each a will be stored in data 
files FOR021.DAT through FOR030.DAT. The Grid and Flow data can be 
separated using a simple fortran code that reads the data and then writes the 
data to two separate files. Splgf.f inputs the desired file ID and outputs a Grid 
file and Flow file. 
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a. Sample Input to Ns2.f 



MACH, 


ALFAO, 


ALFAl, 


ALFARE, 


REDFRE, REYNOLDS 


0.300, 


0.00, 


15.54, 


0.00, 


0.01272, 2.70 


ED2X, 


ED2Y, 


ED4X, 


ED4Y, 


ED 


0.00, 


0.00, 


0.030, 


0.030, 


2.0 


DT, 


COUR, 


NITER, 


NEWTIT 




0.0020, 


2100, 


50, 


1 




RSTRT, 


OSCIL, 


RAMP, 


NPER, 


TSHIFT 


true. 


false. 


true. 


3000, 


-0.5 


TIMEAC, 


IMPLBC, 


EXPLBC, 


CIRCOR 




1, 


false. 


true. 


false 




Vise, 


BLTM, 


JKTM, 


RNGTM 




true. 


true. 


false. 


false 




ITEL, 


ITEU 








31, 


131 








lAl, 


IA2, IA3 


. IA4, 


IA5, IA6, 


IA7, IA8, IA9, lAlO 



1250,1350,1450,1550,1570,1590,1610,1630,1630, 1805 
false 



UNSTST (If true Time = 0.,Set TRUE only when unsteady motion starts from steady 
steady state, TRUE for steady-state ok, but for unsteady restarts must be FALSE 
for proper recording of unsteady motion) 



Mach 

AlfaO 

Alfal 

Alfare 

Redfre 

Reynolds 



Free stream Mach number 

Angle of attack, also mean angle of attack for unsteady 
Amplitude of Oscillatory motion 
Reference angle 

Reduced frequency k = omege * c / 2U 
Reynolds number Re = cU/n 



ED2x 


; X-direction 


2nd 


order 


explicit 


smoothing 


( 


e2x 


= 


0.00 


subsonic. 












0.25 


< 


e2x 


< 


0.50 


transonic) 


ED2z 


: Z-direction 


2nd 


order 


explicit 


smoothing 


( 


e2z 


= 


0.00 


subsonic. 












0.10 


< 


e2z 


< 


0.20 


transonic) 


ED4x 


: X-direction 


4nd 


order 


explicit 


smoothing 


( 


e4x 


— 


0.03 


subsonic. 
















e4x 




0.05 


transonic) 


ED4z 


: Z-direction 


4nd 


order 


explicit 


smoothing 


( 


e4z 


= 


0.03 


subsonic. 



ED 

ISPEC 



Scaling of Implicit smoothing 
Spectral radious parameter 



Dt 

Cour 

Niter 

Newtit 



Time step 

Cour ant number Cu = dt * L_rnax 

Number of Iteration in this run 

Newton siibiteration within each timestep 



RSTRT 

OSCIL 

RAMP 

NPER 

TSHIFT 



Restart 

Oscillatory motion A(t) = AO + A1 * sin { k * M t ) 

Ramp motion 

Number of time steps in one period of oscillation, dt:=T/Nper 
Time shift in radians to start oscilation for any a(t) 



TIMEAC 

IMPLBC 

EXPLBC 



Time accurate Tacc=l and for Jacobian Scaled Dt, Tacc=0 
Implicit wall BC Treatment 
Explicit wall BC Treatment 
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Vise 

TURBL 

JKTM 

RNGTM 


: Only laminar viscosity 

: Baldwin-Lomax eddy viscosity 

: Johnson-King eddy viscosity 

: RNG eddy viscosity 


I TEL, 
lAl etc, 


ITEU : Lower and Upper trailing edge I locations 

: Write out (IAn/100) angles of attack in units 
61 - 70 



Read grid from unit fort. 11 and the flow from unit fort. 31 
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b. Sample Output from Nsl.f 



Grid dimensions : 161x64 



Mach = 


0.30000 










Ao 


0.00000 A1 = 


16.00000 


k = 0.01272 






Re = 


0.2700E+07 Dt = 


0.01829 


Cu= 15.00000 






Dimpl = 


2.00000 










D2x = 


0.00000 D2z = 


0.00000 








D4x = 


0.03000 D4z = 


0.03000 








Restart = 


= T 










Oscil 


= F 










Ramp 


= T 










Itel = 


31 Iteu = 131 


He = 81 








Timeac 


= T 










BC impl 


= F 










BC expl 


= T 










Cour = 


4289.735274387 










L max = 


234539.9275225 










Dt 


0.01829000 










It 


drmax 


dumax 


dvmax demax 


ir 


kr 


1950 


0.849290E-06 0. 


107934E-05 


0.190543E-05 0.174200E-05 


119 


63 


Ni= 1 


0.113997E-06 0. 


302620E-06 


0.448250E-06 0.196797E-06 






it = 1950 


Time = 35.665500 


omega = 0 , 


,007632 phase = 0.043322 






Alfa(t) = 


15.596 Cl = 1 


.425268 Cd = 


= 0.137214 Cm = 0.003384 






2000 


0.778403E-06 0. 


107038E-05 


0.191827E-05 0.158330E-05 


118 


63 


Ni= 1 


0.113176E-06 0. 


297453E-06 


0.450476E-06 0.197668E-06 






it = 2000 


Time = 36.580000 


omega - 0 . 


.007632 phase = 0.044433 






Alfa(t) = 


15.996 Cl = 1 


.453050 Cd = 


= 0.146989 Cm = 0.001525 
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B. DIRECTIONS FOR THE EXECUTION OF U2DIIF.F 



1. Unsteady Case 

A. A complete description of the input and output parameters can be 
found in Teng. 



TABLE 5. U2DIIF.F INPUTS/OUTPUTS 



FILE 


COMMENTS 


FOR001.DAT 


Input 


Uldiif.out 


Output 



B.) For a Unix based machine the code is executed as follows: 

U2diif > U2diif.out 

Note: The original Fortran code can be modified to store desired output in 
separate data files by inserting a WRITE (9,*) statement inside the required do- 
loop. (Example: Cp vs x/c data could be stored in a data file FOR009.DAT). 
This requires some knowledge of Fortran programming. 
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a. Sample Input to U2diif.f 

4 

***************************************************************K** 

AIRFOIL TYPE : NACA 0012 AIRFOIL 
NLOWER = 50 , NUPPER = 50 

0, 50, 50 
12 

0.00, 15.54, 10.67, 0.00, 0.25, 0.00, 0.00 
0 . 00 , 0 . 00 , 0.00 
11.00, 0.010, 0.0001, 0.00 

108, 294, 487, 584, 802, 891, 1006, 1007, 1283, 1554 
ITITLE 

INPUT TITLE (ITITLE = NO. LINE) 

IFLAG NLOWER NUPPER 
AIRFOIL TYPE 

ALPI DALP TOON FREQ PIVOT UGUST VGUST 
DELHX DELHY PHASE 
TF DTS TOL TADK 

ial ia2 ia3 ia4 ia5 ia6 ia7 ia8 ia9 ialO 
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b. Sample Output to U2diif.f 



DATA FROM fort.l: (IRIS) OR FOR001.DAT: (STARDENT) 

AIRFOIL TYPE : NACA 0012 AIRFOIL 
NLOWER = 50 , NUPPER = 50 



IFLAG (0:NACA, 1: INPUT) 


= 


0 


NO. PANELS UPPER SURFACE 


= 


50 


NO. PANELS LOWER SURFACE 


= 


50 


INITIAL ANGLE OF ATTACK 


= 


0.0000 


FINAL ANGLE OF ATTACK 


= 


15.5400 


RISE TIME 


= 


10.6700 


REDUCED FREQ. FOR OSCIL 


= 


0.0000 


PIVOT POINT 


= 


0.2500 


STREAMWISE GUST VELOCITY 


= 


0.0000 


PERPENDICULAR GUST VELOCITY 


= 


0.0000 


X AMPLITUDE OF TRANS OSCILL. 


= 


0.0000 


Y AMPLITUDE OF TRANS OSCILL. 


= 


0.0000 


PHASE OF TRANS OSCILL. 


= 


0.0000 


FINAL TIME 


= 


11.0000 


INITIAL TIME STEP FOR S.S. 


= 


0.0100 


TOLERANCE FOR CONVERGENCE 


= 


0.0001 


FACTOR BY WHICH DTS ADJUSTED 


= 


0.0000 



COORDINATES OF AIRFOIL NODES 



X/C Y/C 

1.000000 0.000000 
0.999013 -0.000141 
0.996057 -0.000562 



0.984292 0.002222 
0.991144 0.001258 
0.996057 0.000562 
0.999013 0.000141 
1.000000 0.000000 



AIRFOIL PERIMETER LENGTH = 2.039290 

STEADY FLOW SOLUTION AT ALPHA = 0.000000 

J X(J) Q(J) GAMMA CP(J) V(J) 
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c 



1 


0.999507 


-0.112073 


0.000000 


0.433845 


-0.752433 


2 


0.997535 


-0.120842 


0.000000 


0.341896 


-0.811236 


3 


0.993600 


-0.126091 


0.000000 


0.279513 


-0.848815 


4 


0.987718 


-0.129448 


0.000000 


0.232223 


-0.876229 


5 


0.979910 


-0.131624 


0.000000 


0.193468 


-0.898071 


6 


0.970208 


-0.132948 


0.000000 


0.160181 


-0.916416 


7 


0.958651 


-0.133600 


0.000000 


0.130705 


-0.932360 


8 


0.945283 


-0.133699 


0.000000 


0.104047 


-0.946548 


9 


0.930159 


-0.133328 


0.000000 


0.079564 


-0.959393 


10 


0.913336 


-0.132558 


0.000000 


0.056812 


-0.971179 



90 


0.894883 -0.131445 


0.000000 


0.035461 


0.982110 


91 


0.913336 -0.132561 


0.000000 


0.056813 


0.971178 


92 


0.930159 -0.133333 


0.000000 


0.079566 


0.959393 


93 


0.945283 -0.133704 


0.000000 


0.104048 


0.946547 


94 


0.958651 -0.133607 


0.000000 


0.130706 


0.932359 


95 


0.970208 -0.132956 


0.000000 


0.160183 


0. 916415 


96 


0.979910 -0.131634 


0.000000 


0.193471 


0.898070 


97 


0.987718 -0.129461 


0.000000 


0.232226 


0.876227 


98 


0.993600 -0.126109 


0.000000 


0.279518 


0.848812 


99 


0.997535 -0.120869 


0.000000 


0.341902 


0.811232 


100 


0.999507 -0.112063 


0.000000 


0.433845 


0.752433 


0 M 


PARISON OF 


G A M M 


A S 





GAMMA FROM KUTTA CONDITION: 

GAMMA FROM CONTOUR INTEGR (TRAIL EDGE) : 
GAMMA FROM CONTOUR INTEGR (MIDPOINTS) : 
GAMMA FROM BOX INTEGR (OFF THE CONTOUR) : 
GAMMA FROM PRECISE CONTOUR INTEG (6 PT) : 



-0.00000024 

-0.00000025 

-0.00000025 

-0.00000024 

-0.00000024 



CD = 0.000324 CL = -0.000001 CM = 0.000001 

*** BEGIN UNSTEADY FLOW SOLUTION **** 



TIME STEP TK = 0.010000 TK - TKMl = 0.010000 



ALPHA(T) = 0.000041 OMEGA(T) = -0.000143 

U(T) = 0.000000 V(T) = -0.000036 



NITR VXW VYW 

0 1.000000 0.000000 

1 0.834730 0.000244 

2 0.827773 0.000255 

3 0.827452 0.000256 



WAKE 

0.010000 

0.008347 

0.008278 

0.008275 



THETA GAMMA 

0.000000 -0.240067E-06 
0.000292 0.436698E-05 

0.000309 0.412307E-05 

0.000309 0.411209E-05 



CONVERGED SOLUTION OBTAINED AFTER NITR = 3 
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J 


X(J) 


Q(J) 


GAMMA 


CP(J) 


V(J) 


PREV PHI ( J) 


CURR PHI(J, 


1 


0.999507 


“0.107774 


0.000004 


0.432728 


-0.753023 


0.032660 


0.032661 


2 


0.997535 


-0.117297 


0.000004 


0.341094 


-0.811709 


0.033045 


0.033045 


3 


0.993600 


-0.123077 


0.000004 


0.279152 


-0.849194 


0.033663 


0.033661 


4 


0.987718 


-0.126843 


0.000004 


0.232351 


-0.876535 


0.034408 


0.034404 


5 


0.979910 


-0.129345 


0.000004 


0.194095 


-0.898321 


0.035211 


0.035205 


6 


0.970208 


-0.130932 


0.000004 


0.161293 


-0.916622 


0.036018 


0.036011 


7 


0.958651 


-0.131801 


0.000004 


0.132278 


-0.932531 


0.036785 


0.036776 


8 


0.945283 


-0.132080 


0.000004 


0.106051 


-0.946691 


0.037475 


0.037463 


9 


0.930159 


-0.131864 


0.000004 


0.081967 


-0.959514 


0.038055 


0.038041 


10 


0.913336 


-0.131225 


0.000004 


0.059576 


-0.971280 


0.038497 


0.038482 


90 


0.894883 


-0.132662 


0.000004 


0.032373 


0.982025 


0.038778 


0.038794 


91 


0.913336 


-0.133893 


0.000004 


0.054048 


0.971077 


0.038497 


0.038511 


92 


0.930159 


-0.134798 


0.000004 


0.077163 


0.959272 


0.038054 


0.038067 


93 


0.945283 


-0.135322 


0.000004 


0.102044 


0.946404 


0.037474 


0.037486 


94 


0.958651 


-0.135406 


0.000004 


0.129134 


0.932188 


0.036785 


0.036794 


95 


0.970208 


-0.134972 


0.000004 


0.159071 


0.916209 


0.036018 


0.036025 


96 


0.979910 


-0.133914 


0.000004 


0.192844 


0.897820 


0.035210 


0.035216 


97 


0.987718 


-0.132066 


0.000004 


0.232097 


0.875921 


0.034407 


0.034411 


98 


0.993600 


-0.129123 


0.000004 


0.279879 


0.848433 


0.033662 


0.033664 


99 


0.997535 


-0.124414 


0.000004 


0.342703 


0.810760 


0.033045 


0.033044 


100 


0.999507 


-0.116362 


0.000004 


0.434960 


0.751843 


0.032659 


0.032658 


: 0 M 


PARIS 


; 0 N OF 


G A M M 


A S 









GAMMA FROM KUTTA CONDITION: 

GAMMA FROM CONTOUR INTEGR (TRAIL EDGE) : 
GAMMA FROM CONTOUR INTEGR (MIDPOINTS) : 
GAMMA FROM BOX INTEGR (OFF THE CONTOUR) : 
GAMMA FROM PRECISE CONTOUR INTEG (6 PT) : 



0.00000411 

-0.00000165 

-0.00000137 

0.00000462 

0.00000385 



CD = 0.000324 CL = 0.005136 CM = -0.003050 



TRAILING VORTICES DATA 



M X(M) Y(M) CIRC 

1 1.004137 0.000001 -0.000009 

TIME STEP TK = 0.020000 TK - TKMl = 0.010000 



ALPHA(T) = 0.000164 

U(T) = 0.000000 



OMEGA (T) = -0.000285 

V(T) = -0.000071 



NITR VXW 
0 0.827452 



VYW WAKE 

0.000256 0.008275 



THETA 

0.000309 



GAMMA 

0.411158E-05 
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CONVERGED SOLUTION OBTAINED AFTER NITR = 0 



J 


X(J) 


Q(J) 


GAMMA 


CP(J) 


V(J) 


PREV PHI ( J) 


CURR PHI 


1 


0.999507 


-0.104926 


0.000011 


0.432928 


-0.753413 


0.032661 


0.032658 


2 


0.997535 


-0.114799 


0.000011 


0.341275 


-0.812045 


0.033045 


0.033042 



20 


0.669285 


-0.109493 


0.000011 


-0.121525 -1.061102 


0.031441 


0.031419 


21 


0.639427 


-0.105975 


0.000011 


-0.138604 -1.069082 


0.029361 


0.029339 


22 


0.609018 


-0.102029 


0.000011 


-0.155875 -1.077078 


0.027013 


0.026992 


23 


0.578179 


-0.097571 


0.000011 


-0.173377 -1.085108 


0.024398 


0.024378 
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C. DIRECTIONS FOR THE EXECUTION OF INCOMPBL.F 
1. Steady State Case 

A. Execution of Incompbl.f is straight forward. First edit the input files, 
FOR001.DAT and incompbl.dat, and set the desired flow conditions. 



TABLE 6. INCOMPBL.F INPUTS 



INPUT 


COMMENTS 


incompbl.dat 


IWAKE-Rag Wake Calulation 
NXT 

NW-No. of point in wake 
ITREND-No. of Iterations 
ITR(l)-Transition Flag Upper Surface 
ITR(2)-Transition Flag Upper Surface 
ISWPMX 

RL-Reynolds Number 
XCTR(l)-Transition Location Upper 
Surface 
IP-Print flag 


FORGO 1. DAT 


Input AOA, pivot point and Airfoil 
Coordinates, and No. of Panels 



B.) For a unix based machine the code is executed as follows: 

Incompbl <incompbl.dat> incompbl.out. 



TABLE 7. INCOMPBL.F OUPUT 



OUTPUT 


COMMENTS 


incompbl.out 


Flow data 







Note: The original Fortran code can be modified to store desired output in 
separate data files by inserting a WRITE (24,*) statement inside the required do-loop. 
(Example: Cp vs x/c data could be stored in a data file FOR024.DAT). This requires some 
knowledge of Fortran Programming. 
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a. Sample Input to Incompbl.f 
File Named fort.l ofFOR001.DAT 

3 

NACA 0012 AIRFOIL 

ALPI PIVOT 

4.000000 0.250000 

NLOWER NUPPER 

50 50 



X/C 



1.00000 


0.98000 


0.96000 


0.94000 


0.92000 


0.90000 


0.88000 


0.86000 


0.84000 


0.82000 


0.80000 


0.78000 


0.76000 


0.74000 


0.72000 


0.70000 


0.68000 


0.66000 


0.64000 


0.62000 


0.60000 


0.58000 


0.56000 


0.54000 


0.52000 


0.50000 


0.48000 


0.46000 


0.44000 


0.42000 


0.40000 


0.38000 


0.36000 


0.34000 


0.32000 


0.30000 


0.28000 


0.26000 


0.24000 


0.22000 


0.20000 


0.18000 


0.16000 


0.14000 


0.12000 


0.10000 


0.08000 


0.06000 


0.04000 


0.02000 


0.00000 


0.02000 


0.04000 


0.06000 


0.08000 


0.10000 


0.12000 


0.14000 


0.16000 


0.18000 


0.20000 


0.22000 


0.24000 


0.26000 


0.28000 


0.30000 


0.32000 


0.34000 


0.36000 


0.38000 


0.40000 


0.42000 


0.44000 


0.46000 


0.48000 


0.50000 


0.52000 


0.54000 


0.56000 


0.58000 


0.60000 


0.62000 


0.64000 


0.66000 


0.68000 


0.70000 


0.72000 


0.74000 


0.76000 


0.78000 


0.80000 


0.82000 


0.84000 


0.86000 


0.88000 


0.90000 


0.92000 


0.94000 


0.96000 


0.98000 


1.00000 




Y/C 

-0.00126 


-0.00403 


-0.00674 


-0.00938 


-0.01196 


-0.01448 


-0.01694 


-0.01935 


-0.02170 


-0.02399 


-0.02623 


-0.02842 


-0.03056 


-0.03264 


-0.03467 


-0.03664 


-0.03856 


-0.04042 


-0.04222 


-0.04396 


-0.04563 


-0.04723 


-0.04878 


-0.05026 


-0.05165 


-0.05294 


-0.05415 


-0.05530 


-0.05634 


-0.05726 


-0.05803 


-0.05868 


-0.05923 


-0.05966 


-0.05995 


-0.06006 


-0.05997 


-0.05966 


-0.05911 


-0.05838 


-0.05737 


-0.05607 


-0.05444 


-0.05236 


-0.04990 


-0.04683 


-0.04309 


-0.03842 


-0.03231 


-0.02382 


0.00000 


0.02382 


0.03231 


0.03842 


0.04309 


0.04683 


0.04990 


0.05236 


0.05444 


0.05607 


0.05737 


0.05838 


0.05911 


0.05966 


0.05997 


0.06006 


0.05995 


0.05966 


0.05923 


0.05868 


0.05803 


0.05726 


0.05634 


0.05530 


0.05415 


0.05294 


0.05165 


0.05026 


0.04878 


0.04723 


0.04563 


0.04396 


0.04222 


0.04042 


0.03856 


0.03664 


0.03467 


0.03264 


0.03056 


0.02842 


0.02623 


0.02399 


0.02170 


0.01935 


0.01694 


0.01448 


0.01196 


0.00938 


0.00674 


0.00403 


0.00126 





File named incompbl.dat 



IWAKE 


NXT 


NW 


I TREND 




0 


161 


37 


40 




ITR(l) 


ITR(2) 


ISWPMX 


RL 


XCTR(l) 


0 

IP 


0 


1 


3000000.0 


0.10000 
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b. Sample Output from Incompbl.f 



NACA 0012 AIRFOIL 



INPUT DATA FOR INVISCID-FLOW CALCULATIONS 

ALPI= 4.0000 PIVOT= 0.2500 

NLOWER= 50 NUPPER= 50 

COORDINATES OF THE BODY 
X/C 



1.000000 


0.980000 


0.960000 


0.940000 


0.920000 


0.900000 


0.880000 


0.860000 


0.840000 


0.820000 


0.800000 


0.780000 


0.760000 


0.740000 


0.720000 


0.700000 


0.680000 


0.660000 


0.640000 


0.620000 


0.600000 


0.580000 


0.560000 


0.540000 


0.520000 


0.500000 


0.480000 


0.460000 


0.440000 


0.420000 


0.400000 


0.380000 


0.360000 


0.340000 


0.320000 


0.300000 


0.280000 


0.260000 


0.240000 


0.220000 


0.200000 


0.180000 


0.160000 


0.140000 


0.120000 


0.100000 


0.080000 


0.060000 


0.040000 


0.020000 


0.000000 


0.020000 


0.040000 


0.060000 


0.080000 


0.100000 


0.120000 


0.140000 


0.160000 


0.180000 


0.200000 


0.220000 


0.240000 


0.260000 


0.280000 


0.300000 


0.320000 


0.340000 


0.360000 


0.380000 


0.400000 


0.420000 


0.440000 


0.460000 


0.480000 


0.500000 


0.520000 


0.540000 


0.560000 


0.580000 


0.600000 


0.620000 


0.640000 


0.660000 


0.680000 


0.700000 


0.720000 


0.740000 


0.760000 


0.780000 


0.800000 


0.820000 


0.840000 


0.860000 


0.880000 


0.900000 


0.920000 


0.940000 


0.960000 


0.980000 


1.000000 





Y/C 

-0.001260 

-0.016940 

-0.030560 

-0.042220 

-0.051650 

-0.058030 

-0.059970 

-0.054440 

-0.032310 

0.043090 

0.057370 

0.059950 

0.056340 

0.048780 

0.038560 

0.026230 

0.011960 



-0.004030 

-0.019350 

-0.032640 

-0.043960 

-0.052940 

-0.058680 

-0.059660 

-0.052360 

-0.023820 

0.046830 

0.058380 

0.059660 

0.055300 

0.047230 

0.036640 

0.023990 

0.009380 



-0.006740 -0. 
-0.021700 -0. 
-0.034670 -0. 
-0.045630 -0. 
-0.054150 -0. 
-0.059230 -0. 
-0.059110 -0. 
-0.049900 -0. 
0.000000 0 . 
0.049900 0. 
0.059110 0. 
0.059230 0. 
0.054150 0. 
0.045630 0. 
0.034670 0. 
0.021700 0. 
0.006740 0. 



009380 -0.011960 
023990 -0.026230 
036640 -0.038560 
047230 -0.048780 
055300 -0.056340 
059660 -0.059950 
058380 -0.057370 
046830 -0.043090 
0.032310 



023820 

052360 

059660 

058680 

052940 

043960 

032640 

019350 

004030 



0.054440 

0.059970 

0.058030 

0.051650 

0.042220 

0.030560 

0.016940 

0.001260 



-0.014480 

-0.028420 

-0.040420 

-0.050260 

-0.057260 

-0.060060 

-0.056070 

-0.038420 

0.038420 

0.056070 

0.060060 

0.057260 

0.050260 

0.040420 

0.028420 

0.014480 



INVISCID FLOW RESULTS 



PANEL 

1 

2 



XP 

0. 99000E+00 
0. 97000E-^00 



YP 

-0.14950E-02 
-0. 44813E-02 



CP 

0.25519E+00 

0.17805E+00 



174 



3 


0.95000E+00 


-0.74415E-02 


0.13179E+00 


4 


0.93000E+00 


-0.10321E-01 


0. 97803E-01 


5 


0. 91000E+00 


-0.13061E-01 


0.74658E-01 


6 


0.89000E+00 


-0.15646E-01 


0.60117E-01 


7 


0.87000E+00 


-0.18113E-01 


0.49465E-01 


8 


0. 85000E+00 


-0.20502E-01 


0.39706E-01 


9 


0.83000E+00 


-0.22829E-01 


0.30916E-01 


10 


0.81000E+00 


-0.25104E-01 


0.21877E-01 



90 


0.79000E+00 


0.27319E-01 


-0.12424E+00 


91 


0.81000E+00 


0.25096E-01 


-0.10630E+00 


92 


0.83000E+00 


0.22823E-01 


-0.88121E-01 


93 


0.85000E+00 


0.20492E-01 


-0. 69791E-01 


94 


0.87000E+00 


0.18084E-01 


-0.49935E-01 


95 


0.89000E+00 


0.15564E-01 


-0.24583E-01 


96 


0. 91000E+00 


0.12899E-01 


0. 93575E-02 


97 


0.93000E+00 


0.10105E-01 


0.51153E-01 


98 


0. 95000E+00 


0.72366E-02 


0. 99432E-01 


99 


0.97000E+00 


0.43440E-02 


0.15874E+00 


100 


0.99000E+00 


0.14480E-02 


0.25519E+00 






INVISCID WAKE 


RESULTS 


PANEL 


XP 


YP 


CP 


101 


0.10033E+01 


0.26712E-04 


0.34783E+00 


102 


O.lOllOE+01 


0.93657E-04 


0.26481E+00 


103 


0.10209E+01 


0.19595E-03 


0.21845E+00 


104 


0.10340E+01 


0.35320E-03 


0.18376E+00 


105 


0.10510E+01 


0.59334E-03 


0.15510E+00 


106 


0.10732E+01 


0.95766E-03 


0.13035E+00 


107 


0.11023E+01 


0.15069E-02 


0.10857E+00 


108 


0.11402E+01 


0.23302E-02 


0.89285E-01 


109 


0.11897E+01 


0.35565E-02 


0.72269E-01 


110 


0.12543E+01 


0.53713E-02 


0.57416E-01 


111 


0.13387E+01 


0.80381E-02 


0.44666E-01 


112 


0.14489E+01 


0.11926E-01 


0.33948E-01 


113 


0.15928E+01 


0.17545E-01 


0.25162E-01 


114 


0.17806E+01 


0.25591E-01 


0.18168E-01 


115 


0.20259E+01 


0.36999E-01 


0.12767E-01 


116 


0.23462E+01 


0.53013E-01 


0.87308E-02 


117 


0.27643E+01 


0.75273E-01 


0.58084E-02 


118 


0.33103E+01 


0.10592E+00 


0.37624E-02 


119 


0.38079E+01 


0.13479E+00 


0.26992E-02 



CL = 0.47937E+00 



INPUT DATA FOR BOUNDARY-LAYER CALCULATIONS 



IWAKE 


NXT 


NW 


ITREND 




0 


161 


37 


40 




ITR(l) 


ITR(2) 


iswpr^ix 


10**-6*RL 


XCTR(l) 


0 

IP 

1 


0 


1 


3.00 


0.10 
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★ ★★★★★★★★★★★★★ CYCXjE 40 ^ 



SUMMARY OF THE DRAG, LIFT AND PITCHING MOMENT 
COEFFICIENTS WITH THE CYCLE 

CD, CL AND CM ARE EVALUATED FROM THE INTEGRATION OF CP 



CYCLE 


CD 


CL 


CM 


1 


0.001978 


0.001684 


0.427162 


2 


0.003184 


0.001730 


0.427756 


3 


0.002734 


0.001684 


0.427162 


4 


0.002422 


0.001730 


0.427770 


5 


0.002227 


0.001684 


0.427157 



35 


0.001684 


0.427157 


0.004733 


36 


0.001731 


0.427790 


0.004494 


37 


0.001684 


0.427165 


0.004731 


38 


0.001731 


0.427807 


0.004491 


39 


0.001684 


0.427167 


0.004731 


40 


0.001731 


0.427822 


0.004488 



BOUNDARY_LAYER PROPERTIES FOR THE LAST CYCLE 

upper SURFACE 

XCTR= 0.134E+00 



NX 


XC 


XS 


CF 


DLS 


UE 


CP 


76 


0.00294 


0.003411 


0.06394 


0.00005 


0.15840 


0.97491 


77 


0.00140 


0.006562 


0.03670 


0.00005 


0.33362 


0.88870 


78 


0.00056 


0.009021 


0.02384 


0.00005 


0.48286 


0.76685 


79 


0.00018 


0.010828 


0.01910 


0.00005 


0.59624 


0.64450 



158 


0.98712 1.018435 


0.00046 


0.00711 


0.87987 


0.22583 


159 


0.99209 1.023461 


0.00035 


0.00749 


0.87320 


0.23751 


160 


0.99639 1.027802 


0.00025 


0.00785 


0.86751 


0.24742 


161 1.00000 1.031453 

LOWER SURFACE — 

XCTR= 0.677E+00 


0.00019 


0.00818 


0.86283 


0.25553 



NX 


XC XS 


CF 


DLS 


UE 


CP 


87 


0.00540 0.000436 


0.47509 


0.00006 


0.02025 


0.99959 



IT 

3 

2 

2 

2 

5 

4 

5 

5 

IT 

2 
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88 


0.00888 


0.004971 


0.15459 


0.00006 


0.18949 


0.96409 


2 


89 


0.01333 


0.010187 


0.02621 


0.00007 


0.36027 


0.87021 


2 


90 


0.01841 


0.016077 


0.01498 


0.00007 


0.53294 


0.71597 


3 



158 


0.98713 


0.992288 


0.00204 


0.00181 


0.89064 


0.20676 


2 


159 


0.99210 


0.997313 


0.00186 


0.00191 


0.88051 


0.22471 


2 


160 


0.99639 


1.001654 


0.00171 


0.00202 


0.87121 


0.24100 


2 


161 


1.00000 


1.005305 


0.00159 


0.00212 


0.86310 


0.25505 


2 
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